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Abstract We live in an age in which high-performance computing is trans¬ 
forming the way we do science. Previously intractable problems are now be¬ 
coming accessible by means of increasingly realistic numerical simulations. One 
of the most enduring and most challenging of these problems is turbulence. 
Yet, despite these advances, the extreme parameter regimes encountered in 
space physics and astrophysics (as in atmospheric and oceanic physics) still 
preclude direct numerical simulation. Numerical models must take a Large 
Eddy Simulation (LES) approach, explicitly computing only a fraction of the 
active dynamical scales. The success of such an approach hinges on how well 
the model can represent the subgrid-scales (SGS) that are not explicitly re¬ 
solved. In addition to the parameter regime, heliophysical and astrophysical 
applications must also face an equally daunting challenge: magnetism. The 
presence of magnetic fields in a turbulent, electrically conducting fluid flow 
can dramatically alter the coupling between large and small scales, with po¬ 
tentially profound implications for LES/SGS modeling. In this review article, 
we summarize the state of the art in LES modeling of turbulent magneto¬ 
hydrodynamic (MHD) flows. After discussing the nature of MHD turbulence 
and the small-scale processes that give rise to energy dissipation, plasma heat¬ 
ing, and magnetic reconnection, we consider how these processes may best be 
captured within an LES/SGS framework. We then consider several specific 
applications in heliophysics and astrophysics, assessing triumphs, challenges, 
and future directions. 
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1 Introduction 

On May 20-23, 2013 a workshop was held at the National Center for Atmo¬ 
spheric Research (NCAR) in Boulder, Colorado, USA entitled “Large-Eddy 
Simulations (LES) of Magnetohydrodynamic (MHD) Turbulence.” The work¬ 
shop was sponsored by NCAR’s Geophysical Turbulence Program (GTP) and 
involved approximately fifty participants from eight countries. 

This review paper is a product of the GTP workshop, though it is not 
intended as a comprehensive account of the proceedings. Rather, it is intended 
as a summary of the issues addressed and the insights achieved, as well as 
an inspiration and a guide to promote future work on this subject. Though 
the subject of interest, namely LES of MHD turbulence, is ostensibly rather 
specific, it encompasses a number of subtle physical processes and diverse 
applications and it draws on the formidable discipline of efficient numerical 
algorithm development on high-performance computing architectures. 

The fundamental challenge that defines the field of LES is that the range of 
dynamical scales active in many turbulent fluid systems far exceeds the range 
that can be explicitly captured in a computer simulation. Examples include 
the convection zones of stars, planetary atmospheres, astrophysical accretion 
disks, and industrial applications such as gas turbines. The central premise of 
LES is that large scales dominate the turbulent transport and energy budget 
so a numerical simulation that captures those scales explicitly will provide a 
realistic depiction of the flow for all practical purposes, provided that the small 
scales that cannot be resolved are somehow taken into account. Strategies for 
incorporating the small scales include explicit subgrid-scale (SGS) models or 
implicit numerical dissipation schemes. 

The range of validity for LES is illustrated schematically in Fig.[^ Gonsider 
a numerical simulation of a turbulent fluid system in which the turbulent 
energy spectrum peaks at some characteristic wavenumber Due to the 
nature of digital computing, any such simulation can only capture a finite 
range in wavennmber, say from L~^ to A~^. Guided by the central premise 
stated above, the lower bound of this wavenumber range often corresponds 
to the largest scales in the system L > £. Meanwhile, the higher bound in 
wavenumber is determined by the resolution limit A which may correspond 
to a numerical grid spacing or to the effective width of some explicit low-pass 
filtering operation that averages over the smaller scales (@- 

If the resolution limit A is smaller than the viscous, thermal, and magnetic 
dissipation scales, collectively represented here as fdiss, then the simulation 
may be regarded as a direct numerical simulation (DNS). However, as noted 
above, DNS are not possible for most turbulent systems in astrophysics and 
space physics. A much more tractable situation is when the resolution limit 
captures the turbulence scale £ but not the dissipation scales; £diss <^A<^£. 
This is the realm of LES. Here SGS models can exploit the self-similarity 
of the turbulent cascade in the inertial range and diffusive prescriptions are 
often sufficient (although even that may not be true in MHD). Ideally, the 
grid spacing A should also be much less than other scales that lead to large- 
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scale anisotropy, such as the Rossby deformation radius, the Bolgiano scale 
of convection, and the pressure and density scale heights. When this is not 
possible, such sources of anisotropy must be taken into account in any explicit 
SGS model. 

Sometimes the characteristic scale of a turbulent flow component is smaller 
than the effective resolution of the simulation or model £< A. One may then 
model the influence of the unresolved scales on mean, resolved flows but one 
might not refer to this as an LES model. A better terminology might be to call 
this a Reynolds-Averaged Navier-Stokes (RANS) approach with some model 
for the turbulent transport, possibly including non-diffusive as well as diffusive 
components. In the following, such calculations are also referred to as mean- 
field simulations (MFS). 

If properly formulated, the LES approach should converge to the DNS ap¬ 
proach as A goes to zero. This is not necessarily true for RANS. For many 
systems such as homogeneous turbulence, there is a smooth transition from 


RANS to LES as the filter scale is decreased from A > i to A i (Schmidt 


20151. SGS models may include non-diffusive transport that resembles the 


Reynolds stress modeling in a RANS system, blurring the distinction between 
the two approaches. Numerical models of such systems may lie anywhere along 
a continuous spectrum of modeling approaches from RANS to LES to DNS. 
On the RANS end of the spectrum, the reliability of the Reynolds stress model 
is paramount, along with analogous prescriptions for turbulent heat transport 
and, in the MHD case, turbulent magnetic induction. For LES models that 
lie more toward the DNS end of the spectrum the details of the SGS model 
presumably become less important, though the simulation becomes more sen¬ 
sitive to the accuracy of the numerical algorithm. Also, as one moves across 
the spectrum from RANS to LES to DNS the computational cost increases 
greatly, along with the number of degrees of freedom. 

In other systems, the transition from RANS to LES is less straightforward; 
one must beware of the “Terra Incognita” that may lie between (Fig. 1). As 
the LES filter size A approaches £, SGS models that rely on the self-similar 
nature of turbulent cascades may break down. There may be a maximum 
scale Ales < £ above which the filtering procedure becomes ill defined and 
unreliable. On the other hand, RANS models may require a sufficient scale 
separation to make statistical averages meaningful, such that Z\ ^ In Fig. I 
this minimum scale for the validity of RANS is labelled as Ameso) in reference 
to the mesoscale modeling of the Earth’s atmosphere. In between these two 
limits, Ales < A < Ameso, lies the Terra Incognita where turbulence modeling 
and simulation can become much more challenging ( Wyngaard||2004 1. 

Due largely to the industrial and atmospheric applications, LES of hy¬ 
drodynamic turbulence is widespread and relatively mature (Sagaut 20061. 
However, most astrophysical and geophysical flows of interest are electrically 
conducting plasmas in which the magnetic field plays an essential dynami¬ 
cal role. For these flows, models must take magnetism into account either 
through the kinetic theory of plasmas (generally necessary for the smallest 
scales) or through the simplifying equations of MHD (often well justified for 
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Fig. 1 The “Terra Incognita” and the realm of LES (from 


Wyngaard||2004 


presented at 


the workshop by Peter Sullivan). Shown is an idealized power spectrum of some turbulent 
field 0 as a function of wavenumber (referred to elsewhere in the paper as k). The peak 
of the spectrum lies at ~ 1/i where ^ is a characteristic length scale of the turbulence. If 
the grid spacing of the simulation, A is much larger than i, then the turbulence is entirely 
unresolved, defining the so-called mesoscale limit of atmospheric science (A > ^meso)- If 
the turbulence is partially resolved, capturing the peak in the spectrum but not the viscous 
dissipation scale id {id A < i)^ then this is the appropriate scenario for LES. 


large scales). Though LES of MHD turbulence can build upon the large body 
of work in hydrodynamic (HD) turbulence, it poses unique challenges that 
must be addressed specifically. These include small-scale anisotropy, nonlocal 
spectral transfer, and magnetic reconnection. 

In section 2 we discuss some of these unique challenges of MHD turbulence 
and highlight particular features of MHD turbulence that may promote the 
development of reliable SGS models. In section 3 we consider the physics of the 
smallest scales where ideal MHD no longer applies, promoting mechanical and 
magnetic energy dissipation and magnetic reconnection, and we ask how these 
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scales may influence the dynamics of the large scales. We then review current 
SGS modeling approaches for MHD in section 4 and assess the triumphs and 
tribulations of current applications in section 5. We summarize the state of 
the field in section 6 and anticipate where it may be headed in the future. 

Though many of the physical processes and challenges we address have 
implications throughout astrophysics, we will focus primarily on solar and 
space physics in this review. This is done to allow us to achieve some depth in 
the material covered while still maintaining a manageable length. For a more 
comprehensive overview of LES/SGS in astrophysics the reader is referred to 


Schmidt (2015). 


2 MHD Turbulence: Challenges and Building Blocks 

2.1 Anisotropy in incompressible unbounded turbulence, from HD to MHD 


Global anisotropy is an essential feature of MHD flows, particularly in the 
presence of a mean magnetic field. The existence of a unique fixed orientation 
yields breaking isotropy towards axisymmetry, with or without mirror symme¬ 
try. System rotation, buoyancy, and density stratification further contribute to 
global anisotropy and inhomogeneity as in HD turbulence. However, in these 
HD cases, if one considers scales small enough (e.g. much smaller than the 
Rossby radius of deformation for rotation) then the constraints become negli¬ 
gible. For MHD turbulence the situation is exactly the opposite; the flow never 
’’forgets” the existence of the large-scale constraint imposed by the magnetic 
field. Indeed as one goes to smaller and smaller scales the anisotropy increases 


(Tobias et al. 2013). This is a severe constraint that must be respected by 


sub-grid scale models. 

In the absence of a mean magnetic held, Alfvenic MHD turbulence can be 


investigated with a more (Iroshnikov 19631 or less (Kraichnan 19651 isotropized 


model. However, even in this case, the substructure of Alfven wave packets 
at small scales cannot be ignored in the overall structure and dynamics of 
the turbulence. If the governing orientation of the small-scale Alfven packets 
is seen as random, a sophisticated stochastic model, mixing anisotropy and 
intermittency, is needed (one such model was discussed by W. Matthaeus at 
the workshop). Again, this is a formidable challenge facing SGS in MHD. 

The theory for MHD turbulence has been developed over the past few 


decades, and there are many recent reviews summarizing various aspects (Kraici- 


2010 


Brandenburg & Nordlund 


2011 


nan fc Montgomery 1980|| Biskamp |2003[ Zhou et al. 2004[ [Petrosyan et al. 


Tobias et al. 2013). One familiar phe¬ 


nomenology is that of interacting wavepackets. This phenomenology arises 
because nonlinear Alfven waves are exact solutions of the full incompress¬ 
ible MHD equations (see, e.g. Parker|[l979 ). A more precise statement is that 
nonlinear interactions only take place when oppositely-signed Elsasser fields 
Z_(. = V + b and Z_ = v — b overlap in space, this statement being valid 
with or without a mean magnetic field Bg, and even in two-dimensional ge- 
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ometry with Bq = 0 where there is no global propagation direction at all. In 
any case one often encounters the heuristic explanation that interactions only 
take place when oppositely propagating wave-packets interact with each other. 
When coherent propagation can occur, it anisotropically interferes with non¬ 
linearity, and gives rise to anisotropic spectra (Shebalin et al.|[l983 Oughton 


et al.|1994). A related effect, the dynamic alignment of turbulent velocity and 


magnetic fields, also has strong effects on MHD turbulence. Global dynamic 


alignment may occur in some ranges of parameter space ( 

Dobrowolny et al. 

1980 

Ting et al. 

1986 

Stribling & Matthaeus 

1991 

Stawarz et al. 

2012 

) as 


a form of long time turbulent relaxation. However local dynamic alignment 
(Milano et al 2001, Boldyrev 2006 Matthaeus et alj 20081 occurs rapidly in 
turbulence. Other types of local relaxation that reduce or suppress the strength 
of nonlinearities imply formation of local patches of correlation associated with 
Beltrami velocity fields and force-free magnetic fields (Servidio et al. 20081. 
Numerical experiments also seem to indicate that the degree of alignment of 
field and velocity is scale-dependent, with the alignment variation even prop¬ 
agating into the dissipative regime ( Boldyrev]|2006 Mason et al.||2006 ). 

For MHD turbulence with a strong externally supported DC magnetic 
field Bo, it is possible to form a large-scale condensate of energy which in¬ 


fluences the turbulent cascade at all smaller scales (Dmitruk & Matthaeus 


20091. Condensation, whether of this type, or of the inverse cascade type, may 


be associated with generation of low frequency 1// noise and long time cor¬ 
relations and sporadic level changes of energy and other quantities over very 
long times Dmitruk & Matthaeus (20071. All of these dynamical effects may 
influence computed solutions, and should be respected by appropriate sub¬ 
grid scale models. These factors, which present a formidable challenge for SGS 
prescriptions, are discussed in more detail in section 

It is possible to describe the second-order correlation tensors with a min¬ 
imal number of correlators, as scalar or pseudo-scalar spectra, accounting for 
the solenoidal properties of both velocity and vorticity fields. The seminal 
stu dies by |Robertson (19351, Chandrasekhar ( 1950| 1951), Batchelor (19821, 
and Craya| (1958) were completed by Oughton et al. (1997) in the MHD case. 
Developed independently by Gambon’s team, a similar formalism improved the 
decomposition in terms of energy, helicity and especially polarization spectra, 
using the orthonormal bases for solenoidal fields, known as a Craya-Herring 


frame of reference (Herring 1974), with its variant of helical modes (Gam¬ 


bon and Jacquin 1989 Waleffe 1992). This formalism is discussed at length. 


with application to turbulence subjected to rotation, density-stratification and 


uniform shear in the recent monograph by 

Sagaut and Gambon 

to 

o 

o 

00 

, and 

extended to the MHD case by Gambon and 

collaborators ( 

Favier et al. 

2012 


Gambon et al. 2012). In addition to the definition of the basic set of spec¬ 


tra and co-spectra, dynamical equations can be written for the correlators, 
generalizing the Lin equation in isotropic turbulence. 

Unfortunately, very few of these results (for both HD and MHD) have been 
used in recent pseudo-spectral DNS in triple-periodic boxes, even if they could 
reproduce anisotropic homogeneous turbulence, despite the finite-box effects. 
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standard discretization, questionable ergodicity from a single realization, and 
other differences with the theoretical context of homogeneous unbounded tur¬ 
bulence. As a first example, helicity cannot be disentangled from directional 
anisotropy (e.g. angle-dependent, or two-component, energy spectrum) and po¬ 
larization anisotropy in DNS started with a single realization, e.g. with ABC 
artificial helical forcing (Salhi et al. 2014). On the other hand, angle-dependent 
spectra and co-spectra, which are not provided by these recent DNS, are useful 
to quantitatively characterize different anisotropic properties, as the horizon¬ 
tal layering in stably-stratified turbulence, and the opposite trend to generate 
columnar structures in flows dominated by system rotation. Such structures 
are often shown only on snapshots in recent DNS, with very indirect linkage to 
statistical indicators, such as one-component, in terms of wavevector modulus 
k or transverse wavevector kj _, spectra. Might there be some analogy between 
the layering in stably stratified turbulence (which is linked to the kinetic en¬ 
ergy cascade of the toroidal velocity component and angle-dependent spectra) 
and the formation of thin current sheets in MHD, as seen in high-resolution 
DNS? 

In addition, the distinction between the 2D ‘vortex’ modes and ’rapid’ 
inertial modes is dependent on discretization in conventional pseudo-spectral 
DNS for purely rotating turbulence, and the dynamics is affected by finite- 
box effects. Only the use of actual confinement, as with rigid boundaries, 
allows one to identify the 2D mode as a dominant one, whereas it is only a 
marginal limit of inertial wave modes in a very large box, and treated as an 
integrable singularity in wave turbulence theory ( Bellet et al.|2006 ). Extension 
of inertial wave turbulence theory, with coupling to ’actual’ 2D modes, was 
recently achieved in a rotating ’slab’ by Scott (2014). 

An important question, useful for SGS modeling in LES, is the range of 
penetration of anisotropy towards smallest scales (see also (3.2). In the HD 
case, an external effect such as mean shear firstly affects the largest scales. 


generating both energy (production) and anisotropy. As suggested by Corrsin 


(1958), isotropy can be recovered at a typical wavenumber, expressed in terms 
of mean shear rate S and the dissipation rate e: ks = \Jje. Similar thresh¬ 
old wavenumbers were proposed by Ozmidov (1965) for stably-stratified tur¬ 
bulence, replacing S by the Brunt-Vaisala frequency, N, and by Zeman (1994) 
in rotating turbulence, replacing S by system vorticity. Even if these simple 
dimensional considerations are only partly supported by DNS or experiments 
( Lamriben et al|[2011 Delache et al.|[20T4 |, they are not sufficient to close the 
problem and to say that anisotropy can be generally neglected at small scale 
in HD turbulence, in contrast with MHD turbulence. Rotating turbulence and 
stably stratified turbulence are much more subtle, because there is no direct 
production of kinetic energy by the Coriolis force, and no direct production of 
total energy, kinetic -I- potential, by the buoyancy force with stabilizing mean 
density gradient (in contrast with turbulence subjected to mean shear). On 
the other hand, a scale-by-scale analysis of the anisotropy in rotating turbu¬ 
lence, without artificial forcing, reveals that the anisotropy first increases with 
increasing wave-number, so that it can be maximum at the smallest scales 
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if the “Zeman wavenumber” ko = is larger than the viscous cutoff. 

These considerations suggest a refined comparison between inertial wave tur¬ 
bulence theory and weak MHD Alfvenic turbulence, with the latter reviewed 
and updated by Boldyrev in the GTP workshop (see Tobias et al.||2013). 


2.2 Is there a need for including advanced backscatter modeling? 


In HD turbulence, backscatter to larger scales plays energetically a significant 
role, but it is usually not systematically correlated with large-scale properties 
of the flow. On the other hand, at least in helical MHD, backscatter plays a 
dramatic role in that it is responsible for the generation of magnetic energy 
at the largest scale through what is known as the a effect. The a effect plays 
therefore an important role in mean-field simulations (MFS), but is ignored in 
LES. 

The a-effect is linked to the upscale transfer of magnetic helicity, which 
occurs in helical MHD turbulence through local (inverse cascade) or nonlocal 
(a-effect) spectral interactions (Bouquet et al 1976; Seehafer 1996; Branden¬ 
burg 2001; Muller et al. 2012). Spectral transfer of cross helicity (it • E) can 
also couple large and small scales and should be taken into account in SGS 
models, as emphasized in the GTP workshop by Yokoi (2013). 

Indeed, cross helicity is produced in the presence of gravity g and a par¬ 
allel magnetic field B, giving rise to a pseudo-scalar g • B that is odd in the 
magnetic field, just like the cross helicity ( Rudiger et al.|2007 |. In such a case, 
a large-scale magnetic pattern emerges, as can be seen from power spectra and 
images shown in Fig. Whether or not this large-scale pattern is a result of 
some inverse cascading of cross helicity, analogous to the a-effect, remains an 
open question; see [Brandenburg et al. (2014| for details. 

Another potential mechanism that may contribute the the generation of 
large-scale structure in MHD flows such as that shown in Fig. is the sup¬ 
pression of small-scale turbulent pressure by a large-scale magnetic field. This 
is currently gaining attention within the context of mean-field modeling of 
the Reynolds and Maxwell stresses. If this suppression dominates over the di¬ 
rect contribution of the magnetic pressure, as is the case for fully developed 
turbulence, then the net effect will be negative (Kleeorin 1989[ Rogachevskii 
20071. In a strongly stratified layer, this can lead to an instability, which is 


now called the negative effective magnetic pressure instability (NEMPI). It 
has been suggested that NEMPI may play a role in causing the magnetic field 
to form flux concentrations in the upper regions of the solar convection zone, 
where the stratification is strongest (Brandenburg et al. 2012; Kapyla et al. 
2012; Kernel et al. 2013). Again, this effect has been successfully captured in 
MFS, where its predictive capabilities have been instrumental in furthering 
our theoretical understanding. We need to ask whether or not LES should be 
modified to include this effect, or whether LES are naturally able to capture 
this type of physics. For example, the dynamic Smagorinsky model used in 


simulations of stellar convective dynamos by Nelson et al. (2011 2013) pro- 
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a: 

H 



Fig. 2 Left: Normalized energy spectra of Bz from an isothermally stratified, randomly 
forced DNS with gjc^ki = 4 (sound speed Cg and forcing wavenumber k{) at turbulent 
diffusive times ^ 0.2 (blue line), 0.5, 1, 2, 5, 10, and 20 (red line). Right: Magnetic field 
configuration at the upper surface near the end of the simulation. Adapted from |Brand enburg| 
let al.lpoTil l. 


motes the generation of coherent flux structures by nonlinear feedbacks that 
are roughly analogous to those responsible for the NEMPI. 


2.3 Possible consequences of misrepresenting the small scales 

A practical goal of LES is clearly to keep the code stable. This means that 
close to the grid scale the flow must become smooth. In reality, the opposite 
is the case: turbulent diffusion of mean flows, mean magnetic fields, mean 
temperature, and mean passive scalars decreases with scales. This picture is 
quite clear in mean-field theory, where turbulent transport coefficients such 
as magnetic diffusivity, rjt, and a effect are known to become wavenumber- 
dependent. We do not yet know whether this plays an important role in LES, 
but we must ask whether certain discrepancies between LES and astrophysical 
reality can be explained by such shortcomings. Below we discuss one such 
example. 

Realistic global dynamo simulations have revealed that magnetic cycles 


are possible at rotation rates somewhat faster than the Sun (Brown et al. 


20111. Yet, the Sun is known to undergo cycles. Could this be a consequence 


of misrepresenting the small scales in the simulations? 

To approach this question, we need to know what is the governing non- 
dimensional parameter that determines the transition from cyclic to non-cyclic 
dynamos. This is a difficult question, because the mechanism behind the so¬ 
lar dynamo is not conclusively identified. Broadly speaking, there are flux 
transport dynamos where meridional circulation plays an important role de¬ 
termining the cycle time and migration of magnetic activity belts. The other 
candidate is just a effect and differential rotation, giving rise to an a-fl dy¬ 
namo in which meridional circulation is unimportant. In mean-field theory. 
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the relative importance of the 17 effect for both scenarios is determined by the 
non-dimensional quantity Cq = Afl/rjtk^. Here rjt is a measure of the turbu¬ 
lent kinetic energy so Cq may be regarded as equivalent to a Rossby number 
based on the differential rotation. The relative importance of the 17 effect over 
the a effect depends on the ratio oi Ca and a similar parameter Ca = ot/'qtk 
that characterizes the strength of the a effect. Both and the ratio CnjCa 
would be underestimated in an LES in which r]t{k)k and a{k) are too big, 
so this suggests that one would need to compensate for this shortcoming by 
increasing 17 to recover cyclic dynamo action. 

Though this reasoning is generally robust, it is based on kinematic mean- 
field theory so its application to MHD LES must be made with care. For exam¬ 


ple, MHD/LES convection simulations by Brown et al. (20111 and Nelson et al. 


(2013) demonstrated a transition from steady to cycling dynamos with both an 
increase in 17 and a decrease in the SGS component of the turbulent magnetic 
diffusivity, rjt- This appears to be consistent with the mean-field arguments in 
the preceding paragraph. However, when the SGS diffusion was decreased in 
the latter case (Nelson et al. 2013), the kinetic energy of the convection in¬ 
creased and the differential rotation weakened due to Lorentz force feedbacks, 
implying an effective decrease in Cq for the cyclic case. Furthermore, though 
the total magnetic energy was greater in the cyclic, low-dissipation case, the 
magnetic topology was more complex, with less energy in the mean (axisym- 
metric) fields. This implies a relatively inefficient a-effect. Further analysis 
confirmed that despite the relatively weak Ail in the cyclic simulation, the 
primary source of toroidal flux was indeed the f7-effect, implying CnjCa > 1- 


3 Small-Scale Dynamics: Dissipation, Reconnection and Kinetic 
Effects 

LES methods as applied to HD have as a central goal the ability to compute 
the dynamics of the resolved scales more accurately under conditions in which 
the Reynolds numbers are too high for a fully resolved DNS computation. Ad¬ 
ditional goals may be to compute the direct influence of the unresolved scales 
on the resolved scales (backscatter), or, in the so-called energy equation ap¬ 
proach, to track the transport of unresolved turbulence. Specific formulations 
of LES relevant to the MHD case will be discussed in more detail in the follow¬ 
ing section. The complications that MHD introduces into small scale physics 
become even more challenging when MHD models are employed to approxi¬ 
mate the dynamics of a low collisionality (or “kinetic”) plasma. These issues 
carry over into additional challenges for LES/SGS modeling. 

In MHD the dissipation function (for simple resistivity and viscosity) is 
known, as in HD. However the phenomenology of MHD dissipation is more 
complex given that both current density structures and vorticity structures 
are available as sites of enhanced heating. Simulations have shown that this 
leads to a dependence of the ratio of kinetic to magnetic energy dissipation 
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2011, 20141, which is 


on the magnetic Prandtl number (Brandenburg 2009 
generally not reproduced by LES. 

Not only is there an additional channel for dissipation, but the nonlin¬ 
ear transfer of energy between velocity and magnetic field is known to be 


more nonlocal in scale than the transfer of total energy across scale (Verma 
2004 Alexakis et al.|20'05). The possibility that these effects come into play in 


varying proportions for MHD flows in differing parameter regimes is not only 
possible but likely, given for example the well known differences in small-scale 
dynamics in the kinematic dynamo regime, the Alfvenic turbulence regime, 
and the large scale reconnection regime. These differences reflect the degree of 
nonuniversality inherent in MHD behavior, which has been demonstrated in a 
number of recent studies ( Lee at al||2010 Wan et ^|2012 1 

Variability in the nature of the cascade represents a challenge in developing 
LES for MHD since the correct modeling of important sub-grid scale physics 
may be situation-dependent. However the challenge is even deeper in the con¬ 
text of low collisionality astrophysical plasmas, even if MHD represents an 
accurate model at the larger scales. At the smaller scales, comparable to ion 
gyroscales or inertial scales, one expects MHD to break down and give way to 
a more complete dynamical description, which is commonly referred to as the 
kinetic plasma regime. Processes occurring at the kinetic scales may resem¬ 
ble analogous MHD processes, but may also differ significantly in their detail. 
The LES developer in such cases may need to understand carefully whether 
the relevant processes to be incorporated into sub grid scale modeling remain 
MHD-like, or if they are possibly influenced strongly by kinetic physics. 

Two examples of processes that are potentially influenced by kinetic physics 
are dissipation of fluctuation energy and magnetic reconnection. These pro¬ 
cesses may lead for example, to electron/ion heating and nonthermal particle 
acceleration. In many space and astrophysical applications of MHD, from the 
solar wind to black hole accretion disks, these mechanisms can play a crucial 
role for the global dynamics of the system, coupling microscopic and macro¬ 
scopic scales. 

One may also ask how the details of small scale processes might have 
influence on the large-scale dynamics that is the emphasis of LES. As long 
as the focus remains in the MHD range of scales, in the usual way one may 
anticipate that energy transfer across scales will be almost independent of 
scale at high Reynolds number. If an accurate estimation of the energy flux is 
available, it enables closure of the SGS problem before the dissipative scales 
are even encountered. This is a key step in the de Karman & Howarth (1938) 
similarity decay hypothesis, and is a familiar component of most HD LES. 
Even the MHD models can become more elaborate, for example when there 
is a need to include backscatter effects (as discussed above). Furthermore, for 
situations that permit inverse cascad^this additional complexity in modeling 
energy transfer becomes mandatory. However when MHD models are employed 


^ Here we distinguish backscatter from inverse cascade, the latter being back transfer, or 
upscale transfer, driven by an additional ideal conservation law. 
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for long wavelength description of kinetic plasma behavior, it transpires that 
there are additional motivations for study of small scale effects in building 
LES models. These potentially include: 

(i) the requirement of following magnetic topology and connectivity, which 
may be influenced by small-scale processes such as magnetic reconnec¬ 
tion, as well as diffusive effects such as Field Line Random Walk (FLRW); 

(ii) the requirement of computing test particle scattering and/or accelera¬ 
tion, in order to employ the models for study of suprathermal particles, 
heat conduction or energetic particles such as cosmic rays and solar en¬ 
ergetic particles; 

(hi) the requirement of representing dissipation, heating and more complex 
kinetic responses (including in some cases radiative cooling), which may 
be regulated by the LES fields. 


In each of the above problems the large scale MHD fields and the cascade that 
they produce establish conditions at the kinetic microscales, and the physically 
significant process - reconnection, heating, particle acceleration, etc., follows 
as a response. It seems clear that an LES model that would include these 
effects must be more elaborate than one that focuses mainly on energy flux. 

Pursuing a better understanding of the small-scale dynamics in MHD tur¬ 
bulence in the inertial range, and even smaller scale kinetic plasma dynamics in 
a turbulent medium, has become a very active area of research in recent years. 
This effort has been boosted by availability of high resolution 3D MHD codes 
and kinetic plasma codes (fully kinetic, hybrid, and gyrokinetics), as well as 
a wealth of new observational data regarding solar wind fluctuations down to 
the electron gyroradius scale ( Alexandrova et ar||2009 Sahraoui et al.pOOQ l. 
These studies have improved our theoretical understanding of the nature of 
the turbulence cascade and its effects as it progresses from magnetofluid scales, 
to proton and electron kinetic scales. The continuation of these advances is 
expected in the next few years to provide a much improved basis for develop¬ 
ment of SGS models that will enable a new generation of MHD and plasma 
LES models. In the following, we shall attempt to provide a brief overview 
of the current state-of-the-art as well as a discussion of key open questions 
regarding small scale dynamics. 


3.1 Do the Small Scales Matter? 

Before taking a detailed look at the small-scale dynamics that must be present 
in any turbulent MHD flow, we must first address a pressing question; do any 
of these details matter? Recall the central premise of LES introduced in 
Since the large scales generally dominate the turbulent transport and energy 
budget, these are the scales we are most interested in; why should we care 
about the small scales at all? 

There are two answers to this question. First, the small scale dynamics may 
influence the large-scale dynamics, often in ways we do not yet understand. A 
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notable example is global MHD simulations of magnetic cycles in convective 
dynamos. Though remarkable progress has been made in recent years, such 
simulations are still quite sensitive to the nature of the SGS dissipation and the 
spatial resolution Charbonneau (20151. This is perhaps not a surprise, since 
the large-scale helds are intimately linked to the small-scale fields by, among 
other things, the topological constraints associated with magnetic helicity. 
Large-scale dynamos rely on these linkages to generate magnetic energy and 
may thus be particularly sensitive to SGS processes. 

More generally, the magnetic connectivity has the distinction of depending 
on microscopic properties such as reconnection activity, while clearly also hav¬ 
ing an influence on the large scale features of the problem at hand. We now 
turn to the solar wind as another example that demonstrates this. During so¬ 
lar minimum conditions the fast solar wind is believed to emanate from polar 
coronal holes while slow wind emerges from nearby regions outside the coro¬ 
nal holes, and perhaps from reconnection activity in coronal streamers. Stated 
this way it is possible that the boundary between fast and slow wind would be 
sharp, but this is not observed; instead the transition is more gradual (Rap- 


pazzo et al 20121. It has been suggested that this boundary is thickened by 
random component interchange reconnection (Lazarian et al|2012b Rappazzo 


et al|20l'2 ) that causes there to be a band of held lines near the boundary that 


have a hnite probability of connecting across this boundary due to dynami¬ 
cal activity. While high resolution codes can simulate small regions near the 
boundary to demonstrate this phenomenon, in an LES scenario it is doubtful 
that resolved scales would contain sufficient information to characterize this 
process. The resolved held lines would be nominal held lines, and if laminar, 
might maintain a sharp boundary at the coronal hole edges. It would be a 
challenge for a rehned LES/SGS model to incorporate sufficient information 
about the space-time structure of the unresolved huctuations so that a model 
could be developed to represent both spatial randomization, due to held line 
random walk, and temporal randomization, due to potentially numerous un¬ 
resolved reconnection sites. 

It is not difficult to hnd other astrophysical plasma problems that depend 
on small scale, or even kinetic scale processes, while also having a signih- 
cant impact on large-scale features. Examples include small and large-scale 
dynamos (15.31 as well as the relative level of electron, proton and minor ion 
heating in the solar wind or in black hole accretion disks. Here, the small-scale 
physics plays a critical role in determining the overall magnetic topology, radia¬ 
tive signatures, and thermodynamics of the system, with signihcant large-scale 
observable consequences. 

The second answer to the question of “why should we care about the small 
sales at all?” is that the small-scale dynamics can potentially have observable 
consequences that are regulated by the large-scale hows and helds. A notable 
example is particle acceleration in solar hares and interplanetary shocks. Sharp 
gradients in large-scale helds promote small-scale reconnection that often pro¬ 
duces a non-thermal spectrum of high-energy particles. These solar energetic 
particles (SEPs) are an important component of space weather, with poten- 
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tial socio-economic consequences. The small-scale reconnection that produces 
SEPs also dissipates energy (and other global quantities) and reshapes the 
magnetic topology. Thus, it may be necessary in some situations to take parti¬ 
cle acceleration into account when devising high-fidelity SGS models. In such 
cases an energy equation formalism would be desirable in order to compute 
particle diffusion coefficients. 

Yet, there are many HD and MHD applications when a simple dissipa¬ 
tive SGS model will suffice. Here the large-scale dynamics is insensitive to 
the small-scale dynamics, provided that the Reynolds and magnetic Reynolds 
numbers are high enough to resolve coherent structures and capture self-similar 
cascades. A notable example here is solar granulation (see ^5.1| ). In this case, 
one would be satisfied with relatively simple LES models, such as ILES (see 

@- 

In order to assess whether or not a sophisticated SGS model is needed, and 
in order to devise such a model when necessary, one must have a comprehensive 
understanding of the fundamental physical processes that operate at small 
scales, and how they influence large-scale dynamics. This is where we now 
turn. 


3.2 Physics of the small-scale cascade 


Laboratory plasmas provided the first quantitative indication that MHD tur¬ 


bulence is anisotropic relative to the large scale magnetic field direction (Robin- 
son fc Rusbridge|1971 Zweben et al. [1979 1, generating spectral or correlation 
anisotropy with stronger gradients transverse to the magnetic field and weaker 
parallel gradients. Simulations in both 2D and 3D demonstrated the dynam¬ 
ical basis for this effect: propagation of fluctuations along the magnetic field 
interferes with parallel spectral transfer, while perpendicular transfer remains 
unaffected ( Shebalin et al.|[T983 Oughton et al.|[l9M |. Gorrelation anisotropy 
of the same type was found to operate relative to the local magnetic field (Gho 
fc Vishniac||2000 Milano et al||200i] ). 


Spectral anisotropy generates a distribution of excitation in wave vector 
such that average perpendicular wavenumbers are greater than average parallel 
wavevectors, i.e., kj_ > fey, relative to the global field. The degree of anisotropy 
becomes greater at smaller scales, so, for example the anisotropy of V X B 


exceeds that of B (Shebalin et al. 19831. Moreover, local correlation anisotropy 
measured by conditional structure functions ( Gho &: Vishniac|2000 Milano et 
al||2001 1 is greater than global anisotropy. 

Another familiar type of anisotropy that emerges in plasma turbulence at 
MHD scales is polarization (or variance) anisotropy. In this case one finds that 
mean square value of each component of the fluctuations perpendicular to the 
mean magnetic field is larger than the mean square parallel component. This 
condition emerges naturally in Reduced MHD treatments of tokamak plasma 


devices, in which the aspect ratio of the device plays a key role (Kadomtsev 


& Pogutse 1974 Strauss 19761 and the resulting nonlinear dynamics is both 
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transverse and incompressible, and also requires spectral anisotropy with fc_L ;?> 
fc|| as discussed above. Later it was shown that Reduced MHD (RMHD) and 
its transverse fluctuations may be derived by elimination of fast magnetosonic 
and Alfvenic timescales in solutions of the full 3D compressible MHD equations 


19921. 


with a strong mean magnetic field (Montgomery 1982 Zank & Matthaeus 


It is noteworthy that the properties of low frequency, high-fcj,, incompress¬ 
ible fluctuations with transverse polarization, equates in wave vocabulary to 
dominance of the oblique Alfven mode, and suppression of the magnetoa¬ 
coustic modes. This characterization of fully developed incompressible inertial 
range MHD turbulence - consisting primarily of a highly oblique spectrum of 
transverse fluctuations has provided a basis for models of plasma turbulence 
by a number of authors (iMontgomery & Turner|19811 |Higdon|1984l iGoldreich 
and Sridhar|[I^ . 


While there are a number of differences in these formulations, they have 
in common that MHD turbulence gets more and more anisotropic at smaller 
scales. One approach ( Goldrelch and Sridhar|19^ ) introduced the term “crit¬ 
ical balance,” to describe the fate of weakly interacting Alfven waves that pro¬ 
duce perpendicular spectral transfer until nonlinear (perpendicular) eddy and 
linear (parallel) Alfvenic timescales become equal. This establishes a relation¬ 
ship between perpendicular and parallel wave numbers that is characterized 
by fc|| oc . As a consequence, one finds fey <C k_i_ at small scales. The same 
relationship is found in the earlier turbulence theory (]Hlgdon||1984| based on 


quasi-two dimensional or RMHD spectral transfer (Montgomery 1982 She- 


balin et aHl983 ) except that the RMHD turbulence energy is mainly confined 


to the region of wave vector space in which the nonlinear time scale Is less than 
the linear wave timescales. The relationship fc|| oc , is common to both, 
if the wavenumbers are regarded as averages of the energy spectrum in the 
inertial range. In any case the preference for perpendicular spectral transfer 


(Shebalin et al. 19831 appears to be a robust result in MHD turbulence and 


should be considered in SGS modeling when there is a uniform magnetic field 
or a very large scale magnetic field present. 

In addition to the energy spectrum, the Inertial range in MHD turbulence is 
characterized by additional correlations. The velocity and magnetic fields are 
typically correlated in direction with the sense of correlation coherent within 


patch-like regions of real space (Milano et al 2001 Matthaeus et al 20081. A 


complementary idea is that the alignment increases systematically with de¬ 
creasing scale ( ]Boldyre~ 2006 Mason et al. 20061. It is also documented 
that turbulence produces patchy, localized correlation of other kinds in MHD 
(IServidio et al.||2008|), and at least some of these appear to be related to the 


tendency for turbulent relaxation (Ting et al. 1986 Strlbllng & Matthaeus 


19911 to proceed locally in cellular regions, such as flux tubes, as a faster, 


Intermediate step towards global decay and relaxation. The types of correla¬ 
tions produced locally and rapidly In this way include (but are not limited 
to), not only the Alfvenic correlation (velocity and magnetic field), but also 
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the Beltrami correlation (velocity and vorticity) and the force free correlation 
(magnetic field and electric current density). All lead to depression of non¬ 
linearity in the inertial range of scales, as seen in the emergence of Beltrami 
correlation in HD ( Pelz et ^|1985 ). 

It is not entirely clear how or whether these additional correlations should 
be included in LES/SGS modeling of the smaller scale MHD cascade. On the 
one hand, the diversity in possible long-term relaxed states suggests dominance 
of different relaxation processes for different parameter regimes. For example, 
to achieve global dynamic alignment, any excess mechanical or magnetic en¬ 
ergy would need to be dissipated. Similarly, in order to achieve global selective 


decay of energy with constant helicity (Montgomery et al. 1978 Matthaeus 


& Montgomery 19801, also known as Taylor relaxation (Taylor 19741, would 


require that mechanical energy be entirely dissipated while magnetic energy 
remains. Presumably these alternative decay prescriptions place different re¬ 
quirements on the nature of dissipation models. Dynamo action with injected 
mechanical helicity at intermediate scales also places requirements on trans¬ 
fer and dissipation rates of energy, magnetic helicity and kinetic helicity (e.g. 


Brandenburg|200fj [Brandenburg &: Nordlund|2011[ [Brandenburg fc Subrama-| 


nian||200~ ). On the other hand if the processes being modeled are principally 

dependent on the decay rate of energy, it may be possible to define energy 
fluxes with relatively simpler prescriptions, such as by partitioning transfer 
between direct and inverse cascade rates. How these issues will influence im¬ 
proved and accurate LES/SGS models for MHD in the future is a current 
research-level problem that is intimately tied in with prospects for universality 
in MHD turbulence, or perhaps universality within classes of MHD behavior. 

While it will likely be necessary to learn more about MHD and kinetic scale 
cascades to build more complete models, it is noteworthy that considerable 
theoretical progress has been made, including computations, by assembling 
turbulence models that may lie somewhat outside of a strictly-defined LES 
concept. These models typically have concentrated on selected effects that are 
thought to be dominant for the chosen problem. Examples of such models are 
mean field electrodynamics ( Moffatt[|l978 Krause fc Radler[|l980 1 often used 
in dynamo theory, Reynolds averaged MHD models such as those used for solar 


wind modeling (Usmanov et al. 20141 and hybrid models based on multiple 


scale analysis and Reynolds averaging (Yokoi et al. 20081. In any of these 


models, we should note that a turbulent resistivity would have essentially the 
same effect as an “anomalous” resistivity, by which we mean a contribution 
to resistivity due to small scale (and also unresolved) kinetic plasma effects. 
This is also an area that has been well studied (e.g. Biskamp]|2000 ). 

As an example of a non-traditional LES approach, Yokoi et al. (20131 de¬ 


scribe a novel self-consistent mean-field theoretical model of turbulent MHD 
reconnection, highlighting cross-helicity dynamo effects. In this, essentially 
sub-grid, model the effects of small-scale turbulence are represented by two 
additional terms in the Ohm’s law: one proportional to the turbulent energy 
density and describing standard effective turbulent resistivity, and the other, 
new term, proportional to turbulent cross-helicity W = {u.' ■ b') and the large- 
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scale vorticity = V X U. Though this model appears to capture the influence 
of small-scale turbulence on large-scale reconnection (Yokoi et al. 


et al. (20f5) found that it does not perform well for supersonic MHD tur¬ 


bulence, where it fails to reproduce the turbulent electromotive force (EMF) 
obtained from high-resolution ILES (standard eddy diffusion models also fail 
in a similar way). More work is needed to determine its viability in different 
circumstances. Indeed, this applies to all SGS models; to the extent that it is 
feasible, their validity and scope should be evaluated by comparing them to 
high-resolution DNS/ILES (e.g. Crete et al.||2015 Meheut et al.|20f5 | and/or 
to kinetic plasma simulations. 

Having introduced some prominent features of the physics of MHD tur¬ 
bulence at small scales, we will now focus on some recent findings from re¬ 
spective studies that are relevant to the fate of the MHD cascade at smaller 
scales. A central issue for many applications in space and astrophysics is how 
the cascaded energy is actually dissipated, and in some astrophysical systems, 
eventually radiated away. For the present purposes the notion of dissipation 
may be described as the irreversible conversion of large scale or fluid scale 
energy into microscopic kinetic degrees of freedom. Important questions that 
have been recently addressed in this area include the response of test par¬ 
ticles to MHD electromagnetic fields, kinetic effects including dissipation of 
cascaded MHD fluctuations and the response in the form of heating, and the 
role of magnetic reconnection, current sheets and tearing, and the associated 
macroscopic effects of changes in magnetic topology and connectivity. 


3.3 Energization and transport of test particles 


The most primitive model of kinetic response to MHD-scale fields is given by 
the test-particle approximation in which the trajectory of individual plasma 
particles is assumed to be determined by the Newton-Lorentz force law, ne¬ 
glecting all feedback of the particle motion on the rest of the plasma or on 
the electromagnetic fields. The basic physics of acceleration, scattering and 
transport, especially of suprathermal and energetic particle populations, is of¬ 


ten discussed in a first approximation using a test particle approach (e.g. Bell 
1978 Jokipii|[r966 1. Not only are test particle studies useful in understanding 
energy dissipation, but in some cases, e.g., cosmic rays and solar energetic 
particles, it is the response of the test particles to the large scale fields, and 
the subgrid scale fields, that is the essential output of the research. 

A self-consistent model extending beyond test particles is needed for accu¬ 
rate representation of the effects on dissipation of slower populations of plasma 
particles, say, those moving at a few Alfven speeds or less. Nevertheless, in spite 
of its shortcomings, the test particle approach, implemented in concert with 
MHD computations, has been valuable for investigation of potential mech¬ 
anisms of energization and dissipation prior to emergence of computational 
capacities that enable equivalent self-consistent kinetic modeling. 
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A good example of this is the use of test particles in the elucidation of the 
role of reconnection and turbulence in energization of suprathermal particles. 
Spectral methods, having favorable resolution properties for turbulence, were 
able to describe the interplay of test particle energization and nonlinear recon¬ 
nection at a relatively early stage ( Ambrosiano et al.|1988 |. In the presence of 
strong fluctuations, reconnection does not settle in to smooth solutions antic¬ 
ipated from tearing mode theory, and instead remains unsteady and bursty, 
and when the Reynolds number at the scale of the dominant current sheets 
exceeds a few hundred, the fluctuations lead to multiple small magnetic flux 


structures, or secondary islands (Matthaeus & Lamkin 1985, 1986 Biskamp 


19861. This subject is revisited in more detail below in section 3.5 


Here we note simply that such structures can entrain or temporarily trap 
test particles, and are strongly associated with the most efficiently energized 
particles. This entrainment and energization was found to occur between mag¬ 
netic X-points and 0-points, as was later found in much greater detail and 
realism using high resolution kinetic plasma codes (e.g. Drake et al.||2006 ). 

It is clear that even a simple model employing MHD simulation fields and 
test particles can begin to identify kinetic effects beyond simple energization. 
Studies showed that small gyroradius particles (e.g., electrons) tend to be 
accelerated in the direction along the electric current sheets, that is, paral¬ 
lel acceleration, while heavier particles (protons, etc) are energized in their 
perpendicular velocities (Dmitruk et al. 2004). Self-consistent kinetic simu¬ 


lations also were able to find this effect, and in fact it is now understood 
through plasma simulation that the regions in and near current sheets are 
sites of enhanced kinetic effects such as suprathermal particles, temperature 
anisotropies, large heat flux, and in general non-Gaussian features of the pro¬ 
ton distribution function (e.g. Servidio et al.|[20T2 Karimabadi et aL]|2013 |. 

More recent test particle studies that employ weakly 3D RMHD simula¬ 
tions ( Dalena et al.|[2M4 | suggest that energization of a single species of test 
particle progresses through at least two stages in the presence of a strong guide 
field with nearly two dimensional low frequency fluctuations: First, at lower 
energies the particles are energized in their parallel velocities, and mainly while 
entrained near reconnection sites inside of current sheets and in essentially in 
accord with the classical neutral point acceleration mechanism. This is also 
sometimes called “direct acceleration.” As suprathermal energy grows and the 
gyroradii become larger than the typical thickness of the current sheets, the 
energization of test particles continue, but with enhancement of perpendicular 
velocities. This has been described as a “betatron” process associated with 
an inhomogeneous perpendicular electric field found near to, but outside of 
strong reconnection sites (Dalena et al. 2014) . This test particle result pro¬ 


vides more detail on earlier results on acceleration in turbulence (Dmitruk et 


al][2004l p^andran||2010D . 


The marriage of test particle studies with high resolution MHD simulation 
has led to a number of insights and questions that are of relevance to ongoing 
efforts to develop SGS models for MHD and plasmas. For example: 
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• When are the processes of cascade, reconnection, test particle energiza¬ 
tion, and dissipation, related? While it is fairly clear that a broad-band cascade 
requires reconnection to occur at various scales along the way, the fact that test 
particles respond to the associated inhomogeneities suggests that dissipative 
processes may be intertwined in this process. 

• The topology of the magnetic field becomes fuzzy when there are nu¬ 
merous small secondary islands, so trapping, reconnection, coalescence and 
particle energization will in general not be explicitly resolved in an SGS/LES 
scheme for a large system. Most likely these features will need to be understood 
well enough to develop a statistical or phenomenological model. 

• If tracking test particle populations at a statistical transport level remains 
a scientific priority in an LES context, then a requirement will be to follow 
parameters needed for SGS transport models, such as SGS energy, character¬ 
istic length scales, and possibly spectral features such as anisotropies, e.g., to 
capture possible resonances. 


3.4 Kinetic effects, dissipation processes and heating 


Once cascading fluctuations reach scales at which kinetic effects become im¬ 
portant, MHD is no longer applicable. In practical terms, this means that for 
an ion-electron plasma, when the cascade arrives at scales as small as either the 
ion gyroradius scale, pi, or the ion inertial scale di = VaI^cp, kinetic effects 
become important and even dominant. For wavenumber k the corresponding 
kinetic range is indicated by kpi > I or kdi > I. To retain effects like finite 
Larmor radii and Landau damping in this regime, one has to employ a kinetic 
descriptior0 Since in space physics and astrophysics we are often dealing with 
low density plasmas for which the collisionality is very weak, it is important 
to keep in mind that some type of effective “collisions” will inevitably cause 
departures from an idealized model such as the Vlasov-Maxwell equations. 
Fundamental effects such as an increase of system entropy and relaxation to¬ 
wards thermal equilibrium, will rely on the presence of these formally small 
contributions to the kinetic equations ( Klimontovich|19^ Schekochihin et al. 


20091. In any case, given the enormous computational cost of nonlinear kinetic 


simulations in six phase-space dimensions, such studies often fall into the cat¬ 
egory of “extreme computing.” Results on turbulence energy dissipation and 
relaxation in turbulent plasmas, employing particle-in-cell (PIC) Vlasov code 
and Eulerian Vlasov codes, are just starting to appear in the literature (e.g. 


Daughton et al.|[201I[ [Servidio et al.|[20T^ [Karimabadi et al.]|2013[ [Haynes ^ 


al.||2014D . 

In light of the impressive continued growth of supercomputing power as 
we head towards the exascale era, it may be expected that kinetic simulations 
will be at the forefront of research into the fate of cascaded energy at kinetic 
scales in the years to come. Meanwhile, various reduced models are being used 


^ By kinetic description we mean a dynamical description of the plasma that involves only 
the one-particle distribution function which depends on velocity, position, and time. 













LES of MHD Turbulence 


21 



Fig. 3 (Left) Volume rendering of magnitude of current density J in the same small region 
of a high resolution 3D MHD simulation at four different time s, showing complex spatial 
structure and evolution in time (adapted from |Mininni et al|2008^ . (Right) Volume rendering 
of J from a 2048® PIC simulation of plasma turbulence, in a periodic box of side 83.8 thermal 
proton gyroradii. Again, fine scale structure is evident, now at kinetic scales. (Courtesy of 
V. Roytershteyn, to be published). 


to complement fully kinetic studies. These include hybrid (fluid electrons and 
kinetic ions), gyrokinetic and gyrofluid models, and an array of fluid models 
that contain some kinetic effects (Hall MHD, multifluid. Finite Larmor radius 
MHD, etc.) 

Fully kinetic 3D PIC simulations of turbulence and reconnection have re¬ 
cently revealed interesting details about kinetic response to cascading MHD 
scale fluctuations. In 3D, ion-scale current sheets spontaneously develop tur¬ 
bulence through various instabilities, producing a chaotic 3D magnetic field 
structure ( Daughton et al.pOll ). Examples of fine scale current structures in 
high-resolution MHD and kinetic simulations are shown in Fig. Interest¬ 
ingly, however, both the reconnection rate and the mechanism for breaking 
the frozen-flux law seem to be unaffected, being close to those obtained in 2D 
simulations. Numerical experiments with different types of initialization, such 
as velocity-shear-driven kinetic turbulence ( Karimabadi et al.|2013 ) show that 
current sheets and intermittency can form in both 2.5D PIC simulations as 
well as 3D PIC simulations. Furthermore, kinetic activity of various types, in¬ 
cluding heating, appears to be concentrated near sheet-like current structures 
(see also Servidio et al. 2012 Wan et al. 20121, which lends credence to the 
emerging idea that a large fraction of kinetic heating may occur in or near cur¬ 
rent sheets and related structures. (For the MHD analogue of intermittency 
associated with current sheets, see section 3.4 below.) Some wave activity is 
also identifiable in the above examples, although in terms of the partitioning 
of energy, this seems to be the exception rather the rule when broad-band 


turbulence develops at kinetic scales ( 

Parashar et al 

2010 Verscharen et al 

2012 Karimabadi et al. 

2013 

|. Kinetic scale complexity and intermittency is 
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also observed using recent high resolution observations in the solar wind (Perri 
eFalllM^ [Wu et al.||2013D . 


Even with recognition of this complexity at kinetic scales, efforts continue 
to analyze the type of wave mode, or perhaps several types, that might be 
viewed in some sense as the elementary excitations from which the turbulence 
is constructed. There are actually several approaches that may be described 
as a wave turbulence approach. On the one hand there is the formal weak 
turbulence theory (e.g. Galtier et al. 20001 that considers the cases in which 
the leading order dynamics is that of propagating waves that obey, to a first 
approximation, a dispersion relation that assigns a frequency to each wavevec- 
tor. Another view is that the distinguishing characteristics of wave modes 
are the polarizations and wave vector directions, which suffice to establish an 
identification. 

Substantial research has been devoted to models that are built upon the 
premise of wave modes that couple to produce wave turbulence in the nonlinear 
regime. One class of such models, defined in the MHD regime, is the Goldreich- 
Sridhar theory (Goldreich and Sridhar 1995| ), (also called “critical balance” 
after one particular assumption that is made in the theory; see Sec. 3.21. 
This theory assumes that all possible excitations are Alfven modes, having 
polarizations strictly perpendicular to the applied mean magnetic field. The 
usual argument is that other wave modes evolve semi-independently so that the 
evolution of the Alfven waves and their mode-mode couplings can be computed 
independently of the magnetosonic wave turbulence. This idea is also routinely 
carried over to kinetic regimes, in which it is assumed that distinct wave modes, 
such as kinetic Alfven waves (KAW), or whistlers, will evolve independently. 
Numerical evidence is usually invoked to support this assumption, but the 
idea remains somewhat controversial. For example, 2.5D kinetic simulations 
that are initiated with Alfven modes, i.e., zero parallel variance, appear to 
generate parallel variance fairly rapidly, although at a lower level. Thus, in 
wave terminology, magnetosonic mode turbulence is generated by Alfven mode 
turbulence within the time span of the current generation of simulations which 
are relative short due to finite computational resources, typically less than, say, 
1000 proton cyclotron periods. Furthermore it is well known that the parallel 
variance component of solar wind turbulence is small but nonzero, as in the 
famous “5:4:1” observations by Belcher & Davis (1971) using Mariner data. 

Goldreich-Sridhar turbulence, which is purely Alfven mode, evolves from 
a wave state through standard weak turbulence couplings towards the critical 
balance state, provided that the zero frequency modes (purely 2D nonpropa¬ 
gating fluctuations) are absent or very nearly absent. This results in a wave 
turbulence that is highly oblique, with mainly near-perpendicular wave vec¬ 
tors involved in the dynamics. At small scales approaching the kinetic range, 
the oblique Alfven modes in Goldreich-Sridhar theory naturally go over to Ki¬ 
netic Alfven waves ( Hollweg]|1999 ), which have received substantial attention 
recently in solar wind observations ( Bale et aT]|2005 Sahraoui et al.|[2010 | 

Another wave mode discussed in connection with wave turbulence and a 
possible role in the kinetic range of solar wind dynamics is the whistler mode 
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( Hughes et al.|[2M4 ), which generally is at higher frequencies than the KAWs 
and probably has lower amplitude in the solar wind, but may still play a role 
in the operative dissipation mechanisms. 

Much of the debate concerning relative roles of wave modes has taken place 
in the context of recent high resolution measurements of fluctuations in the 


dissipation range of solar wind turbulence (Bale et al. 2005 Alexandrova et 


al. 

2009 

Sahraoui et al. 

2009 

2010 


20101. First, it was observed that the electric 


and magnetic field fluctuations as well as the density fluctuations in the scale 
range of <C fc_i_ <C display power law (not exponential) spectra. While 
the knee at k±pi ~ 1 was originally attributed to some form of damping, 
e.g., proton cyclotron damping or Landau damping of kinetic Alfven waves 
(KAWs), it was later suggested that the observed power law exponents can 


be explained solely on the basis of dispersion effects at these scales (Stawicki 


et al. 2001). Consequently, this scale range is now sometimes also called a 


“dispersion range”. Neglecting cyclotron absorption at the proton resonance, it 
may be that significant ion/electron dissipation sets in, respectively, at sub-ion 
scales, kj_pi ^ 1, and electron scales, /cj_Pe ~ 1. Beyond the electron scales one 
might expect the occurrence of exponential spectra ( Alexandrova et ^|2009 ) 
although there are also observations consistent with yet another power law 
( Sahraoui et al.|[2009 ). Clearly, the physics of this entire sub-ion scale range is 
of interest to understanding the heating of turbulent space and astrophysical 
plasmas. The relative importance of ion and proton kinetic mechanisms that 
might give rise to dissipation is likely determined by turbulence amplitude 
( Wu et aT]|2013 ) in addition to kinetic plasma parameters such as the ion-to- 
electron temperature ratio r and the plasma beta /3. 

Standard approaches to studying the physics of the kinetic range of tur¬ 
bulence are Lagrangian PIC and Eulerian solutions of the Vlasov equation, as 
well as the hybrid (fluid electron) variants of each of these. However, there have 
been special reduced models that have emerged that include interesting sub¬ 


sets of the relevant physics. Nonlinear gyrokinetic theory (see Schekochihin et 
al.|2009 and references therein) has been developed in the context of magnetic 


confinement fusion research since the early 1980s, and today it serves as the 
workhorse for computations in tokamak research. An adaptation of gyrokinet- 


ics, as embodied e.g. in the GENE (Jenko et al. 2000) and AstroGK (Howes 


et al. 2008) codes, includes a subset of possible gyrokinetic effects, and has 


been proposed as a model for turbulence investigations in weakly collisional, 
strongly magnetized space and astrophysical plasmas from the inertial range 
through the ion and electron kinetic ranges. The main limitation of standard 
gyrokinetic theory is that it is based on a low-frequency (compared to the par¬ 
ticles’ cyclotron motion) ordering, assuming a decoupling the fast gyrophase 
dependence from the slow gyrocenter dynamics. Notably, this version of gy- 
rokinetics lacks cyclotron resonance and therefore maintains particle magnetic 
moments, thus placing it at odds with some theories of solar wind and coronal 
heating. We note in passing that this constraint can be removed if necessary, 
leading to extended versions of gyrokinetics ( Qin et al.|2000 ). In any case, the 
formulation does include the physics of kinetic Alfven waves, and goes over 
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Fig. 4 (a) Magnetic and electric field spectra in the solar wind obtained from in situ 


measurements by the Cluster spacecraft, from 

Bale et al.| ||2005[|. Corresponding spectra 

from (b) gyrokinetic simulations (|Howes et al. 

|2008|l, and from (c) electromagnetic PIC 


simulations ( [Karimabadi et al.|2013^ are also shown for comparison. In frames (a) and (b), 
wavenumbers are normalized by the thermal ion gyroradius pi \ in frame (c) the wavenumbers 
are normalized by the electron inertial scale de = cjijjpe = VAji^c-e- 


to Reduced MHD in appropriate limits. As such gyrokinetics is well suited to 
describe the Goldreich-Sridhar cascade. Gyrokinetics provides a computation¬ 
ally efficient method to study certain problems, and it has been argued that it 


does capture the physics needed to describe the observed turbulence (Howes 
et al.|[2008l and heating (TenBarge et al.||2013l; this however remains a topic 


of lively discussion. It is also worth noting that gyrokinetics itself has been 
treated using an LES approach ( Morel et al.||20ir 2012). 

Gyrofluid theory is an attempt to reduce gyrokinetics to a multi-fluid ap¬ 
proach via calculating moments and closing the resulting hierarchy of equa¬ 
tions by providing suitable closure schemes ({Hammett et al.|1990 Passot et al. 


2012). Kinetic effects like finite Larmor radii and linear Landau damping can 


be retained with reasonable accuracy, provided that the closures are carefully 
constructed. Also pioneered in the context of magnetic confinement fusion re¬ 
search in the 1990s, gyrofluid models have more recently been tailored and 
applied to various space and astrophysical problems. The development and 
refinement of this approach is a subject of on-going research. 

It is at present unclear which model or models will provide what is needed 
for development of effective LES for low collisonality plasmas. If we knew 
which processes were important, then selection of the appropriate reduced de¬ 
scription models such as 2.5D kinetic codes, hybrid codes with fluid electrons, 
gyrokinetic codes or gyrofluid codes, may provide the efficiency needed to ar¬ 
rive at the needed answers more quickly. However if those processes need to 
be identified, then more demanding 3D fully kinetic Vlasov or PIG represen¬ 
tations may be required. 

We close the section with a few remarks concerning the prediction of turbu¬ 
lence spectra in the kinetic range and the associated ambiguities of identifying 
dissipation mechanisms. Within a wave turbulence framework, one may estab¬ 
lish a contrast between expectations of KAW-mode turbulence and whistler¬ 
mode turbulence. Solar wind observations indicate a dispersive effect near the 
ion inertial scale that has been associated with KAWs (Bale et aL]|2005 1 and 
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with gyrokinetics (Howes et al. 2008); see Fig.|^ However this can qualitatively 


be explained by any dispersive processes that include a Hall effect (Matthaeus 


et al. 2008), so this is not a conclusive observation. On the other hand, the 


KAW and whistler dispersion relations, examined in detail in the observations, 
seem to favor KAWs more than whistlers in kinetic wave number ranges above 
kpi ^ 1 ( Sahraoui et al.|20'09 ). However a strong conclusion cannot be derived 
from this concerning dissipation processes because, first, we do not know that 
the dissipation occurs at these wavenumbers - if it were to occur at much higher 
electron wavenumbers, there may be whistlers or other higher frequency waves 
that actually do the dissipating; and second, we do not know for certain that a 
wave turbulence treatment is even appropriate. Figure [^compares solar wind 
observations with both gyrokinetic and kinetic simulations. In all cases the 
magnetic spectra break to something steeper than the inertial (—5/3) range 
at or near ion kinetic scales. The latter may be the thermal ion gyroradius pi 
or, especially at low plasma beta, the ion inertial scale di = c/uipi = VaI 
(Note that, when the plasma /3 is unity, pi = di ) The PIC result shows a simi¬ 
lar spectrum (possibly closer in slope to —8/3 than to —7/3) with wavenumber 
normalized by electron inertial scale de = c/wpe- Note that d^/dt = \Jme/rrii. 
Evidently the spectra themselves do not strongly differentiate between models. 
Additional examples from computations are shown in Matthaeus et al (20081, 


Sahraoui et al. (2009), and Alexandrova et al. (2009). 


If stochastic acceleration and scattering of orbits (Dmitruk et al. 2004 


Chandran|20I0 Dalena et al.|20I4 ) in or near current sheets absorbs substan¬ 


tial fluctuation energy, then one may have to look beyond linear damping of 
wave modes for an effective dissipation mechanism. Other heating processes 
are also evidently available near turbulent reconnection sites and other co¬ 
herent structures that are seen in very high resolution kinetic simulations 
(Karimabadi et al. 2013 Wan et 2012). Furthermore, other kinetic insta¬ 
bilities (such as firehose and mirror instabilities) and the fluctuations they 


produce may be particularly active near coherent structures (Servidio et al. 
2014) and sharp gradients (Markovskii et al.|2006), and these may contribute 


to dissipation. More work will have to be done to investigate, confirm and 
refine any of these emerging scenarios, in particular with respect to the role 
of non-Maxwellian velocity space features seen in fully kinetic nonlinear sim¬ 
ulations ( Hellinger fc Travmcek|[20l3 Servidio et al.|[20T2 2014). 


3.5 Magnetic reconnection, current sheets and intermittency 


Magnetic reconnection is an important element of the dynamics of plasmas, 
and has been widely studied for more than half a century both theoretically 
and in applications. The basic theory has been well reviewed in terms of MHD 


theory (Biskamp 2000 

Forbes & Priest 2000 

^weibel & Yamada 

2009 
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kinetic effects begin to become important ( [Yamada et al|[^010 ). We will not 
attempt to reproduce such reviews here. Many of the computational models 
in which reconnection has been studied (see references above) have been for¬ 
mulated to include kinetic effects, but in small domains with simple boundary 
conditions and smooth initial data, often based upon an equilibrium current 
configuration. In such cases the reconnection that is studied may be viewed 
as spontaneous reconnection^ which itself is a rich and well-studied approach. 
However it has become increasingly clear that the traditional way of studying 
magnetic reconnection via standard equilibrium current-sheet setups needs to 
be complemented by investigations of reconnection in more complex environ¬ 
ments, and in particular in self-consistent turbulent environments. As of now, 
we are still in the early stages of understanding this more complex situation, 
in particular with respect to kinetic treatments of quasi-collisionless systems. 
One key question in this context is the degree to which phenomena associated 
with reconnection, such as particle acceleration, can act as an alternative route 
to dissipation. Here we can provide only a brief overview of this active research 
topic, which has been examined from a variety of similar approaches. 

The first description of turbulent reconnection (Matthaeus & Lamkin 1986), 
an initial value problem consisting of a sheet pinch evolving in the presence 
of a specified broad-band spectrum of fluctuations, might properly be called 
“reconnection in the presence of turbulence.” In this case one finds bursty, 
nonsteady reconnection, in which one observes sporadically forming intense 
current sheets and vortex quadrupoles, as well as transient multiple X-points. 
This approach was found to lead to elevated rates of reconnection, for resis¬ 
tive MHD, the increase for large systems being comparable to or greater than 
the increase due to Hall effect (Smith et al 2004). An important dynamical 
feature of this problem is the amplification of the turbulence due to nonlin¬ 


ear instability of the initial configuration and subsequent feedback (Lapenta 


2008). 


Another approach, which has usually been applied in three dimensions, also 
begins with a sheet pinch initial condition, but instead of supplying turbulence 
through an initial spectrum of fluctuations, a random source of fluctuations 
acts continuously through a forcing function applied to the region surrounding 
and including the current sheet (see [Kowal et al 2009 Loureiro et al 2009 


Lazarian et al 2012 and references therein). This too gives rise to strong 


turbulence effects, and as might be expected, the reconnection rate generally 
is tied closely to the imposed turbulence amplitude. 

Still another approach in understanding the relationship between turbu¬ 
lence and reconnection is to initialize a system with a large number of magnetic 
flux tubes, as well as random velocity fields, such that the initial state triggers 
a complex cascade. The turbulent dynamics leads to interactions between var¬ 
ious pairs of adjoining magnetic flux tubes (or magnetic islands), leading to 
reconnection with widely distributed reconnection rates, studied in 2D MHD 
(Servidio et al|2009 2010) and more recently in Reduced MHD (see below and 


Wan et al|20I4 ). This type of “reconnection in turbulence” might be viewed as 


similar to both problems described above, but with the random perturbations 
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caused by the cascade itself, rather than being imposed by an initial spec¬ 
trum or by a prescribed forcing function. Some have argued that for systems 
having many flux tubes, this might be viewed as a more natural way to drive 
reconnection with turbulence, but it has the disadvantage of requiring a large 
system, both to establish a high Reynolds number cascade, and to adequately 
resolve the smaller scale current sheets. 


A further complication in understanding reconnection is that its geome¬ 
try can become quite complex in three dimensions ( Priest fc Schrijver||1999 1, 
departing strongly from the familiar two dimensional forms. Nevertheless it 
seems rather certain that in 3D models, as in 2D models, coherent electric cur¬ 


rent structures, including sometimes complex sheet-like structures (Mininni et 


al 2008), continue to play an important role. For example, current sheets in 


RMHD models occupy a central role in models of coronal heating that have 
been studied ( Einaudi fc Velli]1999 Dmitruk et al|1998 Rappazzo et al|20T0 ) 
in the so-called nanoflare scenario. This model is typically viewed as an im¬ 
plementation of the Parker problem ( Parker|1972 ) in which coronal field lines 
are stirred from below by photospheric motions, which causes a braiding or 
tangling of flux tubes, the formation of current sheets between pairs of them, 
and subsequent bursty reconnection and heating. Recently there has been fur¬ 
ther progress in understanding the local statistics of current sheet dynamics 
in the weakly three dimensional Reduced MHD model discussed above, thus 
advancing out understanding of the role of these current sheets in reconnec¬ 
tion, heating and intermittent dissipation in 3D. Other, simpler, models can 
be constructed that take into account the phenomenology of MHD by stipu¬ 
lating that dissipative structures are current and vorticity sheets and that the 
typical time of energy transfer to small scales is modified in MHD when taking 


into account the role of Alfven waves (see, e.g. Grauer et al. 1994 Politano 
and Pouquet||1995[). 


Zhdankin et al (2013) carried out a quantitative statistical analysis of cur¬ 


rent sheets that emerge in RMHD turbulence, and reported on the distribu¬ 
tion functions of current sheets with respect to their dimensions, peak current 
densities, energy dissipation rates and other characteristics. Wan et al (2014) 
reported a similar study using an RMHD coronal heating model, confirming 
many of the results in Zhdankin et al. (2013), while also computing the dis¬ 
tribution of reconnection rates and demonstrating the statistical connection 
between current sheet dimensions and characteristic turbulence length scales. 
An interesting result obtained in both the above studies is that the locus of 
maximum dissipation rate, always the peak of current density for scalar resis¬ 
tivity, is not always, or even usually, located at the component X-points. This 


is a property also found in laminar asymmetric reconnection (Cassak & Shay 


2007), having different magnetic field strength on the two flow sides of the 


reconnection zone. Perhaps not surprisingly, larger reconnection rates are well 
correlated with the proximity of the current maximum to the X-type critical 
point. These purely spatial analyses were extended into the temporal domain 
by Zhdankin et al (2015) who tracked the evolution of dissipative structures 
over time and measured the statistics of their lifetimes and total energy dis- 
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sipation. The results obtained so far for reconnection in RMHD are clearly 
valuable in providing some information about the 3D case, and are specifically 
applicable to low plasma beta, highly anisotropic systems driven at low fre¬ 
quencies, such as the coronal flux tube problem. However the general 3D case 


will undoubtedly contain significantly greater complexity (see e.g. Mininni et 
^|2008 1, much of which remains to be explored. 


The dynamics of the formation of current sheets and other small scale 
coherent structures is of great importance in understanding the intermittent 
cascade and its fate. Furthermore, current sheet formation may be quite dif¬ 
ferent in a large turbulent system than it is in a laboratory device in which the 
magnetic field to leading order is large scale, laminar, and controlled by exter¬ 
nal coils. For example, it is well known that ideal-MHD flows that develop in 
turbulence give rise to intense thin current-sheet structures (Frisch et al||1983' 


Wan et al 2013). The ideal process of current sheet generation, observed at 
short times in high resolution MHD simulations ( Wan et al||2013 ), apparently 
gives essentially identical higher order magnetic increment statistics as are 
seen in comparable high Reynolds number simulations - so we can understand 
that intermittency and the drivers of the conditions that lead to reconnection 
are ideal processes. In retrospect, this could have been anticipated in Parker’s 
original discussion of coronal flux tube interactions ( Parker|[l972 1. Reconnec¬ 
tion may subsequently be triggered at these sites, resulting in dissipation of 
turbulent magnetic field. 

The process of current sheet formation in the presence of weak dissipa¬ 
tion can also involve multiple magnetic X-points and secondary islands or flux 
tubes. This was originally observed in reconnection with finite background 
turbulence, and in that context it was suggested that secondary islands might 


elevate the reconnection rate (Matthaeus & Lamkin 1985 1986 Loureiro et 


al|2009|. Biskamp (1986) suggested that secondary islands might emerge due 


to a linear instability of thin current sheets above magnetic Reynolds num¬ 
bers of about 10^. More recently this instability was revisited based on the 
recognition that the tearing instability can become much faster if it originates 


in a current sheet already thinned to the Sweet-Parker thickness (Loureiro et 
al|2007 Bhattacharjee et al. 2009). The occurrence of this “plasmoid insta¬ 


bility” has been supported by simulation studies using MHD and PIC codes 


(Samtaney et al. 

2009 

Daughton et al. 

2011 

Loureiro et al 

2012 

)• 

Loureiro 

et al 

(2007 

) found that the growth rate of the secondary tearing instability 


of a Sweet-Parker reconnection layer is higher than the inverse global Alfven 
transit time along the layer. M. Velli and collaborators (e.g. Pucci & Velli 


2014) have argued that this result means that, in reality, such a layer cannot 


be formed in the first place. In fact when the the linear growth rate of tearing 
instability equals the inverse of = LjVa the current sheet necessarily is dis¬ 
rupted. This occurs ai a/L where S is the global Lundquist number, 

S = LVa/t]. In the asymptotic limit S —)■ oo, this value is much greater than 
for the SP layer, which suggests that as a current layer is being formed, it is 
disrupted by secondary tearing well before it reaches the SP stage. A similar 
conclusion was presented by (Uzdensky fc Loureiro |[^014|. 
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As mentioned above, one also finds the occurrence of numerous secondary 
islands in a turbulence context, and this will also be the fate of any multiple is¬ 
land scenario at finite amplitude. High resolution 2D MHD turbulence simula¬ 
tions display a proliferation of magnetic X-points at sufficiently high magnetic 
Reynolds number Rm (Wan et al|2012 1, with the number of X-points and flux 
tubes observed in the simulations, scaling as Rm ■ A cautionary word is that 
secondary islands (or plasmoids) also form due to numerical error, and there 
is a requirement for careful resolution studies ( Wan et ffi||2010 1 to ensure that 
complex multiple plasmoid reconnection is physical and not numerical. It is 
unclear at present what precise relationship exists between plasmoid instability 
and generation of secondary islands turbulence. It is noteworthy that simula¬ 
tions, even laminar cases, usually trigger reconnection with a finite amplitude 
perturbation, and by the time multiple islands are observed there are many 
finite amplitude modes participating. Whether the origin is instability or cas¬ 
cade, it seems certain that at high magnetic Reynolds numbers, reconnection 
zones will be complex, even in 2D, so that the details of the many individual 
reconnection processes will almost certainly not be resolved in LES, but rather 
their aggregate effect will need to be built into a SGS model. 

To close this section we note first the implications of the evolving perspec¬ 
tives on reconnection for an application, say, solar coronal heating. In such a 
complex driven system that is far from equilibrium, one should properly view 
current sheet formation, dissipation, magnetic reconnection, and nanoflares, 
not as independent processes, but rather as outcomes of a nonlinear MHD- 
turbulent cascade in a self-organized solar corona. Building such effects into 
an LES/SGS model will be challenging. 

Finally for completeness we list some of the outstanding questions that 
emerge from this discussion: 


1. Do singular (ideal) structures matter for the dissipative case? 

2. Does the 2D case matter to understand the 3D case? 

3. Are rotational discontinuities a central piece of 3D reconnection? 

4. Does current sheet roll-up play a role? 

5. What role do invariants (magnetic and cross helicity) have in reconnection? 

6. Is the rate of dissipation independent of Reynolds number? 

7. What are the dissipative and/or reconnecting structures, and how to iden¬ 
tify them in a real (natural) system? 

8. What is the role of the magnetic Prandtl number? 

9. Do small-scale kinetic effects that emerge in reconnection alter large-scale 
dynamics? How? 

10. Can Adaptive Mesh Refinement help in a general approach? 


4 LES in MHD 

The discussion in the previous section ((|^ highlights some of the challenges 
in devising reliable SGS models for LES of MHD turbulence. In this section 
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we describe in more detail the practical implementation of LES/SGS model¬ 
ing, focusing on explicit approaches that employ formal filtering operations 
designed to decompose the flow into large and small-scale components. The 
large-scale motion is computed by solving the Altered non-stationary equations 
of MHD while the SGS terms are parameterized and expressed in terms of the 
Altered quantities. Though real plasma flows are likely much more subtle 
current MHD-LES models often assume that the subgrid scales (SGS), also 
referred to as subfilter-scales (SFS), are relatively isotropic, homogeneous, and 
universal. 


LES is a method for simulation of flows with large Reynolds numbers. It is 
generally not valid for low Reynolds number flows since it assumes that there 
is a substantial (order unity) nonlinear transfer to small scales. Furthermore, 
as discussed in the usual assumption of isotropy at small scales may not 
be realized. This may occur in rotating and/or stratified flows if the cut-off 
wavenumber where the filter is applied is in the anisotropic range. And, it 
may occur more generally in turbulent MHD flows where anisotropy may pro¬ 
gressively increase toward smaller scales unless this is mitigated by turbulent 
reconnection processes which may help recover isotropy. 

Initially the Large Eddy Simulation technique was developed for the simu¬ 
lation of HD turbulence of neutral fluids, particularly in the context of atmo¬ 
spheric and engineering applications (Meneveau and Katz||2000 Sagaut|[2006 


Glazunov and Lykossov 2003). This has been extended to MHD turbulence 


by several authors who adapted known HD closures to the MHD case and 


developed new SGS models. ( 

Agullo et al. 

2001 

Vliiller and Carati 

2002a|b 

Yoshizawa|I987 

Zhou and VahalaJ1991 

Knaepen and Moin|2003 

). Our discus- 

sion of the general approach follows that given in 

Chernyshov, Karelsky and 


Petrosyan (2006a|. Our intention is to illustrate how the MHD equations can 


be cast into a a traditional LES framework. We make no attempt at a com¬ 
prehensive survey of existing LES/SGS models. For alternative approaches see 


Miki & Menon 

(2008 

), Crete et al.|(2015|), and 

Schmidt (2015). Of particular 

note is the model presented by Crete et al. 

(20151 in which the SGS stress ten- 


sor involves nonlinear correlations between the resolved (filtered) velocity and 
magnetic field gradients, along with Smagorinsky-like expressions for the SGS 
kinetic and magnetic energy proportional to the corresponding rate of strain 
tensors. This is based on a Taylor expansion of the velocity and magnetic field 


within a filter box as originally proposed for HD by Woodward et al (2006a). 


As stated above, a LES applies a filtration operation to the primitive equa- 

For the incompressible MHD equations. 


tions as suggested by Leonard (1974 


the filter G satisfies the following normalization property: 



Xj, Aj)dx'j = 1. 


( 1 ) 


^ For an innovative filtering approach based on wavelet transforms and adaptive mesh 
refinement see |De Stefano &; Vasilyev| | |201^ 
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Here G{xj — Xj,Aj) is the filter itself, of width Aj. Then, the filtered 
velocity is expressed as follows: 


Uj= u{xj)G{xj-Xj,Aj)dxj, ( 2 ) 

J a 

where a = Xj — ^Aj and b = Xj + ^Aj, Aj = (Aa;,Ay,A^). Xj = {x,y,z) 
are axes of Cartesian coordinate system. The other physical fields are filtered 
similarly. 

Let us present all the variables of the problem as the sum of a filtered 
(large scale) and unfiltered (small scale) component: u = u + u',B — B + B', 
p = p + p' etc., with Uj, Bj the velocity and magnetic induction components, 
and p the pressure. 

To simplify the modeled equations describing compressible MHD flows, it 
is convenient to use mass-weighted filtering (also known as Favre filtering) 
so as to avoid the appearance of additional SGS terms. It is determined as 
follows: _ 

/=^, (3) 

P 

with p the density. The over line in eq. ([^, denotes ordinary filtering while the 
tilde denotes mass-weighted filtering. Mass-weighted filtering is used for all 
physical variables other than the pressure, density, and magnetic fields. 

The Favre filtered velocity takes the following form 


puj ^ jl^pUjG{xj-Xj,Aj)dxj 
P la Pi^j)G{Xj - Xj,Aj)dXj 


( 4 ) 


The Favre filtered quantities can be presented in the form of a sum, for instance 
for the velocity: u = u + u", where the double prime designates the small-scale 
component. 


The Favre-filtered MHD equations take the following form (Chernyshov, 
Karelsky and Petrosyan|[2006a I: 


filtered continuity equation 


^ df^ ^ 
dt dxj 

— filtered momentum conservation equation 


dpui 


d 

dxj 


(^pUiUj + pSij 


1 , 52 

Re 


2M2 


— filtered induction equation 


( 5 ) 


dxj ’ 

( 6 ) 
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as 


— where 


rjBj — rjEj = 0, 


^ij ^ij 0 , 




'kij — ^fkSkk^ij 


(SkkS.. 


ij, 


— and p - density; p - pressure; Uj - velocity in direction xj; 
aij = 2pSij — ^pSkk^ij + C^kk^ij - viscous stress tensor; 

S^J = 1/2 {dui/dxj + duj/dxi) - strain rate tensor; 
p - dynamic viscosity; C - bulk viscosity; 

5ij - the Kronecker delta; Sijk - the Levi-Civita symbol; 

rj = IAi:a - magnetic diffusion; a - specific electric conductivity; 

Fi = SijkjjBk/c - Lorentz force; B - magnetic field; j - current density. 

Re = poUqLq/po is the Reynolds number, Rm = uoLq/tjo - the mag¬ 
netic Reynolds number. = uq/Cs the Mach number, Mq = uojua - the 
magnetic Mach number, and Cg is sound speed determined by the relation: 
Cs = V IPo / Po ; Ua is the Alfven speed, Ua = Bq /The bulk coefficient 
of viscosity ^ is neglected. 

The terms on the the right-hand-side of equations §-0 designate the 
influence of the SGS terms on the filtered component: 


t “ = p {u,Uj - UiUj) - {B^Bj - B,Bj) ; 


( 8 ) 


= {uiBj - UiBj) - {BiUj - BiUj) 


(9) 


Note that compressibility alters the form of the SGS stress tensor r/J but 
the magnetic SGS tensor rk is the same as for incompressible MHD. The SGS- 
scale hydrodynamic pressure in typically neglected in the filtered equations for 


compressible, neutral flows with low Mach numbers (Piomelli 1999 Zang et 


al. 1992). By extension, the SGS magnetic pressure is also neglected in eqn. 
(8). However, attempts are being made at more general models. For example. 


Grete et al. (2015) incorporate the full SGS strain rate tensors, including the 


symmetric components that account for SGS kinetic and magnetic pressure. 

Let us consider the filtered equations @-0 in more detail. As far as 
the small-scale velocity (and the other flow variables) u" = u — u \s unknown; 
it has to be estimated with the use of the large-scale velocity obtained by 
means of filtration. Theoretically, there is no functional dependence between 
the small-scale velocity u" and the large-scale one u, so any estimation of u” 
will contain error. DNS can sometimes be used to estimate this error but this 
can only be carried out for relatively low Reynolds numbers due to limited 
computational resources. 
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Thus, the filtered system of MHD equations contains the unknown tur¬ 
bulent tensors r” and The task of the SGS model is to express these 
unknown tensors in terms of the filtered flow components Ui and Bi using 
some sort of turbulent closure (parameterizations). Ideally, the closure model 
should capture such effects as the Richardson turbulent cascade. 

Let us consider closures for r,"- and To guarantee the non-negativity of 
subgrid energy, these tensors must satisfy some conditions, called realizability 
conditions. A necessary and sufficient condition of non-negativity is provided 
by the positiveness of the semidefinite form for the turbulent tensors Tij such 
that: 


Tii >0 for i G {1,2,3} , 

I Tij I < y/TuTjj for i,jG {1,2,3}, 


( 10 ) 


det(Tij) > 0 . 

Let us assume that the form of the turbulent tensor r“ is analogous to the 

. . 


viscous stress tensor (eddy viscosity model), while is analogous to ohmic 


dissipation. This yields: 


1 


1 


T^- - oT^k^ij = -2vt S,j - -SkkSi 


n'^kk^ij ~ 


( 11 ) 

( 12 ) 


where 


Sii = X 
^ 2 


1 f dui du 


dxi 


_ ■£ 

dxi 


is the large-scale strain rate tensor; 

1 fdB, dB, 


J^3 - 2 


dxj 


dxi 


is a large-scale magnetic rotation tensor. Here vt and ? 7 t are scalar turbulent 
diffusion coefficients that may in general depend on the spatial coordinates 
and time. 

In the right-hand side of equations 0 and ( |12[ ) the symmetric terms of 
the magnetic rate-of-strain tensor: 

= [dB.jdxj -b dBjIdx,) /2, 

and vorticity tensor: 

J“ = {dui/dxj — duj/dxi) /2 

are omitted, because their contribution is negligible in many circumstances 
(Muller and Carati 2002a). However, Crete et al. (2015) have incorporated 
these symmetric terms and find that they may be important in particular 
for supersonic flow regimes, as encountered, for example, in the interstellar 
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medium. We note again that the main purpose of SGS modeling is not to 
fully reconstruct the information lost due to filtration but rather to accurately 
capture the influence of the SGS flow on the large-scale energy distribution 
and transport. 

Often the term \T^}^5ij is combined with the thermodynamic pressure, 
V(p -I- where k = (rn -|- T 22 + T33)/2 is the SGS turbulent kinetic 

energy (Erlebacher et al. |T9^ . In the present paper we consider the isotropic 
component explicitly, though the isotropic component of the magnetic tensor 
(121 vanishes because of vanishing Ja). The isotropic component of r“ can be 


found from the realizability conditions (10), which give 


■^12 + '''13 + "^23 ^ T11T22 + T11T33 -I- T22T33 


(13) 


By using we obtain the following expression for the isotropic com¬ 

ponent of 


k>lV3{ut\S-\) 


(14) 


where |^“| = 


1/2 


The anisotropic and isotropic components of 
can then be obtained from Eqs. © and ( [T4| ). 

Different closures for the compressible MHD equations were developed 


by Chernyshov, Karelsky and Petrosyan (2006b | and further analyzed by 


Chernyshov, Karelsky and Petrosyan (2007). LES of MHD turbulence were 


compared with DNS and it was shown that the five closure models considered 
provide sufficient dissipation of kinetic and magnetic energy at comparatively 
low computational expense. 

The effects of heat conduction were considered by [Chernyshov, Karelsky] 


and Petrosyan (2006c 2008a) who developed models for the SGS terms in 


the energy equation as well as for the magnetic terms in the momentum and 
induction equations. LES of decaying MHD turbulence were performed and 
their greater efficiency compared to DNS was demonstrated. SGS models were 
similarly validated for studies of self-similar regimes in forced turbulence by 


Chernyshov, Karelsky and Petrosyan (2010 2012). Further work on the LES 


of compressible MHD turbulence focused on the local interstellar medium 


Chernyshov, Karelsky and Petrosyan 

2008b) and the kurtosis and flatness 

of the turbulent flow Chernyshov, Karelsky and Petrosyan ( 

2009 

). For a com- 

prehensive review of SGS modeling see 

Chernyshov, Karelsky and Petrosyan 


(2014). 


Other models can be developed, using for example expressions deriving 
from the integro-differential equations for energy spectral density in the frame¬ 
work of two-point closures of turbulence like the Eddy Damped Quasi Normal 
Markovian closure ( [Chollet and Lesieur||1981 ). Newer versions have seen sev¬ 
eral developments implemented, such as (see Baerenzung et al. (2011) and 


references therein): (i) including both eddy diffusivities (viscosity and resis¬ 
tivity) as well as eddy noise; indeed, the effect of the small scales on the large 
scales is potentially a dissipation of energy (although eddy diffusivities can 
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be negative) as well as a stochastic forcing; (ii) adapting to spectra that dif¬ 
fer from the classical Kolmogorov law, a feature that can be useful in the 
presence of magnetic fields, rotation or stratification, i.e. when the resulting 
energy spectra may be in a weak turbulence regime due to wave-eddy interac¬ 
tions; and (iii) including in these two types of coefficients the effect of helicity 
(velocity-vorticity correlations) as encountered in tropical cyclones or in the 
Planetary Boundary Layer. Helicity can be created by a combination of ro¬ 
tation and stratification and can play a role in the generation of large-scale 
magnetic fields. 


In fact, in MHD, there are two other helical fields that can be defined, 
namely the cross-correlation between the velocity and the magnetic induction, 
and the magnetic helicity (correlation between vector potential and magnetic 
induction in three space dimensions); their effect on the large-scale dynamics 
of MHD flows has been considered in [Yokoi (2013) where the role of cross¬ 


correlation on turbulent reconnection is particularly stressed (see also Yokoi 
and Yoshizawa|[l9M). 


There are other types of LES that have been developed. Of particular note 
are implicit LES (ILES) methods that do not include any explicit SGS model 
but do include intrinsic dissipation and dispersion due to the nature of the 
numerical algorithm. These methods received much attention at the workshop 
and make up a growing fraction of astrophysical and geophysical turbulence 
simulation. Their primary attraction is maximal resolution; dissipation oper¬ 
ates only at the grid scale, leaving larger scales essentially free of artificial 
diffusion. However, it can be difficult to assess the influence of the dissipation 
scheme on the properties of the resolved flow, particularly in the case of MHD 
where nonlinear spectral interactions are intrinsically nonlocal. For further 
details on ILES see Grinstein et al. (2011) and Schmidt (2015). 

Another possibility that has been tested in two dimensions in MHD in a 
pseudo-spectral code is to decimate Fourier modes after a given cut-off scale 
( Meneguzzi et al.]|1996 ); the rationale behind such a decimation, whereby for 
example half the modes are taken for Fourier shells beyond the chosen cut-off 
(and the method can be iterated), is that there are 47rfc^ modes of characteristic 
wavenumber k, i.e. a large number at small scales; it is to be expected that their 
role is statistical (and with a stochastic component) and that therefore these 
modes do not have to be all treated explicitly. For example, in a computation 
with 3072^ grid points and with a de-aliasing using a 2/3 rule, the ratio of 
maximum to minimum wave numbers is 1024. Of the 27“'' x 10® modes in such 
a computation, half of them (or roughly 13 billions) are for wave numbers 
k > 710, with « 12 X 10® in the very last Fourier shell (of unit width) alone 
( Marino et al.|[2013 ). 

In conclusion, it has been suggested that a possible future methodology 
to tackle complex turbulent flows as found in astrophysics and space physics 
might be to combine multiple approaches, including ILES, explicit SGS mod¬ 
eling, and also adaptive mesh refinement ( Woodward et al||2006a De Stefano 
& Vasilyev |Schmidt|[M5t 
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5 Applications 

In this section we briefly consider several applications of LES in MHD of 
relevance to astrophysics and space physics in order to highlight both the 
successes and the challenges of the field. 


5.1 Realistic LES of solar granulation 


A prime example demonstrating the success of LES is solar granulation. This 
includes in particular the uppermost visible surface of the solar convection 
zone where radiation transport and time-dependent ionization are important. 


1976 

1 on compressible stellar convection in the 1970s, 

Nordlund 

(1982 

1985 


pioneered the field of realistic convection simulations of the solar surface, which 


has since advanced considerably 1 

Stein & Nordlundll 19891 

1998| 

Steffen et al. 

1989 Wedemeyer et al. 

2004 

Vogler et al. 

2005 Rempel 

2014 

). These sim- 


ulations employ a tabulated equation of state together with fully nonlocal 
radiation transfer and realistic opacities. They used different combinations of 
subgrid scale modeling including Smagorinsky viscosity, shock-capturing vis¬ 
cosities, hyperviscosities, Riemann solvers, monotonicity schemes, etc., which 
can be classified as implicit LES (ILES). 

Simulations of solar surface convection reproduce solar observations re¬ 
markably well, both qualitatively and quantitatively. An example is shown in 
Fig.i The intensity contrast in the simulated granules is about 16%, which 
agrees with the observed contrast of about 10% after taking atmospheric see¬ 


ing and the telescope point spread function into account (Stein & Nordlund 


19981. Other quantitative successes include power spectra, spectral line for¬ 


mation, acoustic mode excitation, and local dynamo action (Nordlund et al. 
200^|Rempel|[M4l ). 


An important question is to what extent the success of these simulations is 
due to the ILES technique employed, or to aspects of the physics that make this 
application particularly amenable to LES modeling. For example, the strong 
density stratification makes convection highly anisotropic, with a dilution of 
vorticity tending to “laminarize” upflows while the dynamics of the downflows 
are controlled mainly by buoyancy and entrainment. So, details such as the 
forward transfer of kinetic energy into the dissipation scale are not specifically 
tested. Furthermore, near the visible surface, the radiative diffusivity is large 
enough that no SGS model is needed for the internal energy equation. Thus, a 
key component of the dynamics is effectively captured through DNS. This may 
also account for the success of geodynamo models, which are able to run with 
a realistic value for the magnetic diffusivity, relegating SGS to the velocity 
field alone ( Glatzmaier||2002 Jones|[20TT ). 
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Fig. 5 Comparison between a granulation pattern from a simulation with 12 km grid size 
(left), an observed granulation pattern from the Swedish 1-meter Solar Telescope at disk 
center (middle), and the simulated one after convolving with the theoretical point spread 
function of a 1 meter telescope. The simulation images are for wavelength integrated light 
intensity while the observed image is for a wavelength band in the near UV. The image 
was taken on 23 May 2010 at 12:42 GMT with image restoration by use of the multi-frame 
blind de-convolution technique with multiple objects and phase diversity ( |van Noort et ah 
2005|l. Court esy of V. M. J. Henriques and G. B. Scharmer and adapted from [Brandenbu^ 
&: Nordlund] l [^011| . 


5.2 The bottleneck effect in HD turbulence 


Incompressible forced turbulence simulations have been carried out at resolu- 


tions up to 4096^ meshpoints ( 

Kaneda et al. 

2003 

|. A surprising result from 

this work is a strong bottleneck effect ( 

Falkovich 

1994 

I near the dissipative 


subrange, and possibly a strong inertial range correction of about k~^'^ to 
the usual inertial range spectrum, Interestingly, similarly strong iner¬ 

tial range corrections have also been seen in simulations using a Smagorinsky 
subgrid scale model ( Haugen and Brandenburg||2QQ6 ) using 512^ meshpoints; 
see the dashed line in Fig. Here we also show the results of simulations 
with hyperviscosity, i.e. the diffusion operator has been replaced by a 
Z/ 3 V® operator (Haugen and Brandenburg 20041, also with 512^ meshpoints 
(dash-dotted line). Hyperviscosity greatly exaggerates the bottleneck effect, 
but it does not seem to affect the inertial range significantly; see Fig. 
Woodward et al ( 2006a|b ) found a similar bottleneck effect in ILES of ho¬ 


mogeneous, decaying, sonic turbulence (Mach number 1). Furthermore, they 
implemented a nonlinear SGS model (mentioned previously in designed 
to supplement the numerical dissipation and they found that the combined 
ILES+SGS model could both alleviate the bottleneck effect and reproduce 
the spectrum of higher-resolution ILES across much of the resolved dynamical 
range. 


If the details of the inertial range spectrum are sensitive to the dissipation 
even in this most fundamental of applications then what should we expect for 
more complex flows? Does this call into question the central premise discussed 
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0.001 0.010 0.100 


Fig. 6 Comparison of energy spectra of the 4096^ meshpoints run i Kaneda et al.||2003 i, 
solid line, and 512^ meshpoints runs with hyperviscosity (dash-dotted line) and Smagorin- 
sky viscosity (dashed line). (In the hyperviscous simulation we use 1/3 = 5 x The 

Taylor microscale Reynolds number of the Kaneda simulation is 1201, while the hypervis¬ 
cous simulation of |Haugen and Brandenburg] ( |2004| | has an approximate Taylor microscale 
Reynolds number of 340 < Re;^ < 730. For t he Smagorinsky simulation the v alue of Re>^ is 
slightly smaller. Courtesy of Nils E. Haugen l |Haugen and Brandenburg|200^ . 


in ^ that the dynamics of the large scales can be reliably captured despite 
the challenges of modeling the SGS physics? 


5.3 Problems in dynamo theory 
5.3.1 Small-scale dynamos 

Unlike many industrial applications where LES have been tested against exper¬ 
iments, this is currently impossible for hydromagnetic flows exhibiting dynamo 
action. Except for astrophysical dynamos, there are not even a hand-full of lab¬ 
oratory experiments to date that produce self-excited hydromagnetic dynamo 
action. Therefore, an important benchmark is provided by DNS. 

Dynamos come in two flavors: small-scale and large-scale dynamos (see 
Brandenburg fc Subramanian| [20051 [Brandenburg et al.||2012] for recent re¬ 
views). In the kinematic phase, during which the magnetic field grows expo¬ 
nentially from a weak seed magnetic field, the magnetic field exhibits a 

spectrum. Evidently, this spectrum diverges toward small 
one cannot expect to obtain the correct growth rate with 
LES. This has consequences for understanding the excitation conditions of 
small-scale dynamos. DNS have demonstrated that the onset of small-scale 
dynamo action depends on the value of the magnetic Prandtl number (see 


Kazantsev (1950 


length scales, so 
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Iskakov et al. 2007 and Brandenburg 2011 for the classical incompressible case 
and Fedderrath et al. 2014 for the case of supersonic turbulence). Meanwhile, 
the magnetic Prandtl number does not enter into traditional LES, so this 
question cannot be addressed. This is also the case for the magneto-rotational 
instability (MRI), which has been shown to be sensitive to the value of the 
magnetic Prandtl number. 

As the small-scale dynamo saturates, the peak of magnetic energy moves 
gradually toward larger scales, so there is a chance that this can be modeled 
with LES. However, simulations using a Smagorinsky-like magnetic diffusivity 
prescription have yielded saturation field strengths that are significantly below 
those obtained with DNS ( Haugen and Brandenburg|2006 1. This shortcoming 
might be related to not yet being in the asymptotic regime in which both 
magnetic and kinetic Reynolds numbers are large enough. A similar situation 
might apply to the ratio of kinetic to magnetic energy dissipation, which is 
known to scale with magnetic Prandtl number to a power that is around 
0.6. Again, LES are not currently able to shed any light on this, because the 
magnetic Prandtl number does not enter in standard subgrid scale models. 
Addressing this deficiency is a difficult but perhaps auspicious challenge in 
need of further research. 


5.3.2 Large-scale dynamos 


Large-scale dynamos produce magnetic fields whose scale exceed that of the 
turbulent eddies. They are believed to be relevant for understanding the global 
22-year cycle of the Sun’s magnetic field, and similar large-scale magnetic helds 
in other astrophysical bodies. A leading theory for understanding large-scale 
dynamos is mean-field theory which explains the occurrence of correlation 
in the mean electromotive force that has a component parallel to the mean 
magnetic field. This is generally referred to as the a effect, is typically related 
to the presence of helicity in the system. However, DNS have shown that its 
magnitude is reduced with increasing values of the microphysical magnetic 
Reynolds number (Cattaneo & Hughes 19961. This is now generally referred 
to as catastrophic quenching and has to do with a magnetic contribution to the 
a effect, which is proportional to the current helicity of the fluctuating field, 
j-b. Here, j = Vxb/ is the current density of the fluctuating magnetic field, 
b. It is in turn related to the magnetic helicity of the small-scale field, a ■ b, 
where a is the magnetic vector potential of 6 = V X a. It obeys the evolution 
equation 

d ___ _ 

—a ■ b = —uxb ■ B — 2?7/roJ • 6 — V • (exa), (15) 


Here, uxb is the mean electromotive force, rj is the microphysical magnetic 
diffusivity, and e X a is the magnetic helicity flux, where e is the fluctuating 
electric field. We shall now discuss two aspects of this equation that are relevant 
to LES. 

First, if the system is homogeneous, i.e., there are no boundaries and no 
large-scale variations of turbulent intensity and helicity across the system, the 
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divergence of the magnetic helicity flux vanishes. In that case, there must be 


a balance between the first two terms on the right-hand side of Equation (151. 
Although the assumption of homogeneity is only of academic interest, it does 
provide a test case that must be obeyed equally by DNS and LES. In partic¬ 
ular, replacing the microphysical diffusion by hyperdiffusion changes the scale 
dependence of the j ■ b term and has been shown to exaggerate the amplitude 


of the large-scale field relative to that of the small-scale field (Brandenburg 
fc Sarson]|2002 )■ This phenomenon is well understood, but it would still be of 
interest to experiment with other representations of small-scale magnetic dis¬ 
sipation to see how the large-scale magnetic field is being artificially modified 
by the numerical representation of the small-scale physics. 

Secondly, in the inhomogeneous case, magnetic helicity fluxes are possible. 
In principle, these fluxes are gauge-dependent, and would thus be unphysical. 
If, in the statistically steady state, the time derivative of a • 6 on the left-hand 


side of Equation (151 vanishes, we have 


0 = —2uxb ■ B — 2r]^Qj • b — V • (eXa) 


(16) 


which implies that V • (exa) is now balanced by terms that are manifestly 
gauge-independent. This is a remarkable property that allows us to measure 
the magnetic helicity flux divergence. This has been done in several recent pa¬ 


pers ( Mitra et al.|20ld Hubbard &: Brandenburg|2010 Del Sordo et al.|2013 |. 
Interestingly, it turns out that the V- (exa) remains subdominant for all sim¬ 
ulations performed so far, and that only at the largest resolution available so 
far, it becomes approximately equal to the 2r]^Qj ■ b term. Here, the magnetic 
Reynolds number based on the wavenumber of the energy-carrying eddies is 
about 1000, which is still barely achievable in DNS. 

It is at present unclear whether LES are able to lead to meaningful insight 
into the regime of larger magnetic Reynolds numbers. One practical difficulty 
in determining e X a is the computation of the magnetic vector potential, which 
is not always readily available. 


5.4 Astrophysical Turbulence and MRI 


Most astrophysical plasmas are highly compressible. Thus, capturing the in¬ 
fluence of shocks is essential. Since the viscous scale in shocks is far too 
small to be resolved, DNS is not possible so modelers must turn to LES. All 
shock-capturing astrophysical codes include some sort of subgrid-scale model. 


whether it be an explicit artificial viscosity (e.g. Stone & Norman 1992) or 


an implicit numerical dissipation that arises through a Riemann solver as in 


Godunov-type methods (e.g. Roe 1986 Leveque 20021. These methods have 
been used for over 50 years, and there is no question this approach works 
extremely well. 

Developing more sophisticated SGS models for turbulence in highly com¬ 
pressible astrophysical plasmas will be a formidable challenge. A complete 
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understanding of how shocks and contact discontinuities interact with the tur¬ 
bulent cascade and affect small scales is still lacking and is a major research 
problem in its own right. 

One very important example of astrophysical turbulence in which com¬ 
pressibility does not play an essential role is the magnetorotational instability, 
or MRI ( Balbus &: Hawley|1991 ). MRI-induced turbulence is thought to dom¬ 
inate the energy and angular momentum transport in magnetized accretion 
disks ( Balbus fc Hawle}^|1998 ). 

The MRI is a particularly important case study for LES of MHD turbu¬ 
lence because the role of artificial dissipation has been closely scrutinized in 
recent years. This scrutiny began with an influential paper by [Fromang &:| 


Papaloizou (2007) that demonstrated a decrease in the amplitude of the MRI- 


induced turbulent stresses with increasing spatial resolution for the specific 
case of a local shearing box with no net flux through the layer. The poten¬ 
tial implications were profound; If this trend were to continue to the dynamic 
ranges active in actual accretion disks, then the turbulent transport would be 
drastically less efficient than previously thought and insufficient to account for 
the outward angular momentum transport necessary to sustain the accretion 
process ( Balbus fc Hawley||1998 ). They attributed this behavior to the form of 
the numerical diffusion and its interaction with the source terms that sustain 
the instability. Other LES simulations soon confirmed this result for similar 
model configurations and demonstrated that the decrease in turbulent stresses 
with increasing resolution does not occur when explicit diffusion is included 
(see Bodo et al.||2011 and references therein). 

However, after nearly eight years of active research, it appears that this 
“convergence problem” is not as serious as initially thought; the problem ap¬ 
pears to be a symptom of this particular model setup and goes away when 
other model configurations are considered. For example, if the vertical extent 
of the computational domain is increased so that it is twice as large as the 
horizontal extent, the convergence problem goes away, the turbulent stresses 
increase by an order of magnitude, and the Fourier power spectrum is altered 
substantially. Furthermore, the convergence problem does not arise for shear¬ 
ing boxes with a background density stratification and/or a net magnetic flux 
through the layer (e.g. Froman^|2013[ Turner et al.||2014 ). 

In the no-net-flux case, the MRI must sustain a shear dynamo and it may 
be that the properties of this dynamo (which is a macroscopic flow related to 
the outer scale of the turbulence) is not properly captured in small boxes and 
this introduces an artificial dependence on SGS diffusion. Further insight into 
why the MRI behaves so differently in small, non-stratified, shearing boxes 
with no net flux requires a deeper understanding of the MRI-driven dynamo, 
possibly based on reduced models. This deeper understanding is also needed 


to account for the magnetic cycles found in the larger boxes (see Fromang 
2013 Turner et al.||2014 ). 


Furthermore, direct measurement of turbulent resistivity in the MRI by 
three different groups all find the same result, using very different codes and 
Reynolds numbers (Guan fc Gammie|2009 Lesur fc Longaretti|2009 Fromang 
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fc Stone|2009|. The turbulent magnetic Prandtl number is close to one, which 


implies the resistivity is very large (of order ui, where u and I are characteristic 
turbulent velocity and length scales; this implies that the turbulent magnetic 
Reynolds number is much smaller than that given by the Ohmic resistivity). 
This argues that macroscopic (turbulent) effects are more important than 
microscopic diffusivities, which in turn argues that the MRI is not sensitive to 
SGS physics. 

Further confidence in the ILES approach to modeling MRI turbulence 
comes from the recent study by Meheut et al. (20151. They compared lower- 
resolution ILES to high-resolution DNS and found good agreement, at least for 
the special case of a non-zero mean field and a low magnetic Prandtl number, 
meaning that the magnetic diffusion was captured explicitly while the kinetic 
energy dissipation was relegated to the numerical diffusion. In particular, the 
kinetic and magnetic power spectra at low to intermediate wavenumbers in a 
DNS with resolution (800, 1600, 832) were well reproduced by ILES runs with 
resolutions of order 128^ and 256^ at a fraction of the computational expense. 

If, on the other hand, one wishes to construct explicit SGS models for LES 
of MRI, it may prove beneficial to exploit the potential magnetic induction 
introduced by Salhi et al. (2012), B-V0 where Q is the potential tempera¬ 
ture. This is an analogue of the potential vorticity which, unlike the potential 
vorticity, is a Lagrangian invariant for a magnetized, Boussinesq fluid. 


5.5 Hybrid Kinetic-MHD Models 


In 1 3.4 we emphasized that the small-scale dynamics of a plasma flow are 
often not well represented by MHD. This is particularly the case for low- 
density plasmas such as the solar wind or Earth’s magnetosphere. Departures 
from MHD must be treated by solving the kinetic equations in some form, 
often with simplifying assumptions designed to mitigate the computational 
requirements. See §3.4| for further details and for a survey of some applications 
from solar and space physics. 

Many other astrophysical flows require a kinetic-MHD description to cap¬ 
ture the essential physics. Here there is no question that SGS physics is im¬ 
portant. For example, anisotropic conduction and viscosity on small scales can 
influence the large-scale dynamics through the magneto-thermal and magneto- 
viscous instabilities ( ]Balbus|2000 2001 Quataert et al.|[2002 ). There is some 
hope that hybrid PIC simulations of small scales will be able to provide a rea¬ 
sonable SGS model (including coefficients of conduction and viscosity) that can 
be used in kinetic MHD models of various astrophysical problems including the 
MRI (see j |5.4[ ), hot gas in galaxy clusters, and turbulence in the interstellar 
medium. 

During the GTP workshop, W. Schmidt described a promising hierarchi¬ 
cal approach for including small-scale kinetic reconnection effects as an SGS 
model in MHD simulations. The approach is similar to self-consistent mean- 
field dynamo theory but encompasses three stages: 
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1 . kinetic simulations of reconnection and dissipation providing effective trans¬ 
port coefficients (e.g., the turbulent EMF a and fi coefficients) to account 
for kinetic plasma processes; 

2. In turn, these coefficients are to be used in non-ideal MHD simulations of 
turbulent dynamos and reconnection, providing a sub-grid model for stage 

3; 

3. quasi-ideal simulations of MHD turbulence. 


6 Summary and Outlook 


The diverse applications surveyed in section all have at least two things 
in common. First, these systems are well described on large scales by the 
equations of MHD and second, they are characterized by turbulent parameter 
regimes that are inaccessible to DNS. Computers simply are not capable of 
modeling all relevant scales from the macroscopic scales to the Larmor radius. 
Thus, in order to model such systems it is absolutely essential that we adopt 
an LES approach. 

For many systems it may be sufficient to simply minimize the artificial dis¬ 
sipation through the use of numerical methods that include their own intrinsic 
dissipation. Such methods, often referred to generally as implicit LES (though 

for a more precise definition of ILFS), maximize 


Grinstein et al. 2011 


see 

a simulation’s dynamic range for a given spatial resolution by confining the 
artificial dissipation to scales comparable to the grid spacing. 

Other applications may be more subtle. For these, it may be necessary to 
model the subgrid-scale physics more reliably in order to accurately capture 
the dynamics of the larger scales. This can be achieved by applying a formal 
filtering procedure to the governing equations (Sec. and then introducing 
parameterized or tabulated SGS models based on theoretical and phenomeno¬ 
logical arguments or on local simulations that capture the small-scale plasma 
physics. 

When devising SGS models for MHD, some guidance can be provided by 
the much more mature field of LES/SGS modeling in turbulent HD (non¬ 
magnetic) flows, which has received much attention particularly in the context 
of atmospheric and engineering applications (Sagaut 2006), with a growing 


body of literature also on highly compressible astrophysical flows (Schmidt 


2015). Still, as discussed in sections 2|l3| MHD possesses its own unique chal¬ 


lenges. 

Some of the challenges in representing SGS physics arise from the nature 
of the MHD equations themselves. The presence of magnetism in a turbulent, 
electrically conducting fluid introduces an intrinsic anisotropy to the flow that 
becomes more pronounced with decreasing scale ()|^. The small-scale flow is 
also intrinsically inhomogeneous, marked by intermittent patches of enhanced 
dissipation and magnetic reconnection in current sheets. This small-scale dissi¬ 
pation and reconnection can heat the plasma and reshape the large-scale mag¬ 
netic topology. Other factors that can influence the coupling between large and 
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small scales include magnetic helicity, cross helicity, dynamic alignment, and 
the suppression of small-scale turbulence by large-scale magnetic flux {{2.2). 
Most of these effects are neglected by current SGS models, which often assume 
some degree of isotropy and homogeneity to make paramphenomenaeteriza- 
tions more tractable. Further investigation of these issues and associated SGS 
model development is sorely needed. Help in dealing with inhomogeneity and 
intermittency may also come from adaptive mesh refinement (AMR), which 


was auspiciously addressed at the GTP Workshop by O. Vasilyev (De Stefano 
& Vasilyev |[M^ . 


Yet, the challenges in modeling SGS physics do not stop there. In most 
plasma flows (particularly those with low density), the MHD equations cease 
to be valid on the smallest scales, giving rise to kinetic effects that lie outside 
the scope of ideal or even resistive MHD. Such effects may regulate the dis¬ 
sipation of energy and magnetic helicity, the associated plasma heating, and 
the restructuring of the magnetic topology through magnetic reconnection. 
Kinetic effects also introduce new phenomena that may influence large-scale 
dynamics, including non-thermal particle acceleration and anisotropic heat 
conduction and viscosity. A promising path forward is to couple MHD models 
to kinetic or hybrid codes that capture some of the relevant kinetic effects 


( 13.4 5.5). But this will not be easy; it is a formidable theoretical and com¬ 


putational challenge. 

As mentioned several times in this review, a promising way to validate 
LES models in MHD and to guide their development is to compare them 
with higher-resolution DNS or ILES. Preliminary results from astrophysical 


application have generally been promising (Grete et al. 2015 Meheut et al. 


20151 but more work is needed. Kinetic and hybrid simulations can also be 


used to motivate and assess SGS models, particularly for models that are not 
purely dissipative. We expect to see much progress on these fronts in the next 
5-10 years. 

There is no oracle or omen to tell you whether your application requires 
sophisticated SGS modeling or if ILES is sufficient. This judgement must be 
made on a case-by-case basis grounded on a thorough understanding of the 
underlying physics and indeed, it is still being assessed even for the relatively 
well-established problems surveyed in m Though there are robust features 
of MHD turbulence that can be exploited in SGS models, many aspects of 
the SGS physics are likely not universal. Yet, if prudence is followed when 
designing numerical models and interpreting the results, current applications 
do give us confidence in the central premise of LES ((Q, namely that real 
heliophysical and astrophysical systems can be meaningfully modeled when 
only a fraction of the dynamically active scales are explicitly resolved. LES 
of MHD turbulence is a still a nascent field, brimming with challenges...and 
opportunities. 
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